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ABSTRACT 

We use the halo model of clustering to compute two- and three-point correlation func- 
tions for weak lensing, and apply them in a new statistical technique to measure properties of 
massive halos. We present analytical results on the eight shear three-point correlation func- 
tions constructed using combination of the two shear components at each vertex of a triangle. 
We compare the amplitude and configuration dependence of the functions with ray-tracing 
simulations and find excellent agreement for different scales and models. These results are 
promising, since shear statistics are easier to measure than the convergence. In addition, the 
symmetry properties of the shear three-point functions provide a new and precise way of 
disentangling the lensing i?-mode from the i?-mode due to possible systematic errors. 

We develop an approach based on correlation functions to measure the properties of 
galaxy-group and cluster halos from lensing surveys. Shear correlations on small scales arise 
from the lensing matter within halos of mass M ^ Thus, the measurement of two- 

and three-point correlations can be used to extract information on halo density profiles, pri- 
marily the inner slope and halo concentration. We demonstrate the feasibility of such an anal- 
ysis for forthcoming surveys. We include covariances in the correlation functions due to sam- 
ple variance and intrinsic ellipticity noise to show that 10% accuracy on profile parameters 
is achievable with surveys like the CFHT Legacy survey, and significantly better with future 
surveys. Our statistical approach is complementary to the standard approach of identifying in- 
dividual objects in survey data and measuring their properties. It can be extended to perform 
a combined analysis of the large-scale, perturbative regime down to small, sub-arcminute 
scales, to obtain consistent measurements of cosmological parameters and halo properties. 
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1 INTRODUCTION 

The gravitational lensing of distant galaxy images due to large-scale structure, known as cosmic shear, has been well established as a 
cosmological probe (e.g, Mellier 1999, Bartelmann & Schneider 2001 and Wittman 2002 for reviews). Many independent groups have 
reported significant detections of two-point shear correlations, providing constraints on the mass density of the universe (r^mo) and the mass 
power spectrum amplitude (erg) (e.g., Van Waerbeke et al. 2001b; Bacon et al. 2002; Refregier, Rhodes & Groth 2002; Hoekstra et al. 2002; 
Brown et al. 2002; Hamana et al. 2002; Jarvis et al. 2003). Non-linear gravitational clustering substantially enhances the cosmic shear signal 
on angular scales ^ 10' (Jain & Seljak 1997). An accurate description of non-linear effects on shear correlations is important to interpret 
the measurements. 

The most effective way of studying non-linear structure formation has been TV-body simulations, which are accurate on scales larger than 
the numerical resolution limit, since the relevant physics is only gravity on scales of interest. However, current survey data require accurate 
models of large-scale structure over a huge dynamic range of length scales. A simulation needs to sample cosmological scale (~ lOOMpc) 
in order to have a fair sample. On the other hand, lensing statistics at relevant angular scales are affected by highly non-linear structures, dark 
matter halos with a size ;S Mpc, and structure that needs to be resolved down to scales approaching O.OlMpc. The required dynamic range is 
prohibitive for current computational resources. In addition, to perform multiple evaluations in model parameter space requires thousands of 
simulations runs, which is also prohibitive. While some numerical short-cuts can be used to get around these constraints, having an accurate 
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analytic model would be a valuable complementary, and sometimes essential, tool. It could also provide physical insights into the complex 
non-linear phenomena involved in gravitational clustering. 

The dark matter halo model of clustering fulfills such a niche. The method was first constructed in real space to understand the correlation 
functions of the galaxy distribution in terms of halos (Neyman & Scott 1952; Peebles 1974; McClelland & Silk 1977; Scherrer & Bertschinger 
1991; Sheth & Jain 1997; Yano & Gouda 1999; Ma & Fry 2000a,b; Sheth et al. 2001; Berlind & Weinberg 2002; Takada & Jain 2003b, 
hereafter TJ03b). The model was also formulated in Fourier space, since this leads to simpler expressions for the Fourier-transformed 
counterparts of the higher-order moments (Seljak 2000; Peacock & Smith 2000; Ma & Fry 2000c; Scoccimarro et al. 2001; Cooray & Hu 
2001a,b; Takada & Jain 2002 hereafter TJ02; also see Cooray & Sheth 2002 for a review). Given the model ingredients: the halo density 
profile, mass function and bias, each of which has been well investigated in the literature, the halo model can be used to compute the statistics 
of cosmic fields. There are several advantages in using the halo model. First, it allows us to easily implement multiple model predictions 
in parameter space. Second, the measured signals are explicitly understood in terms of the halo model ingredients. This will be relevant for 
the comparison with other observations such as the X-ray or the Sunyaev-Zel'dovich (SZ) effect for clusters of galaxies. Third, the model 
accuracy can be easily refined by incorporating results from Af-body simulations with higher resolution. Interestingly, so far the halo model 
has led to consistent predictions to interpret observational results of galaxy clustering as well as reproduce simulation results (e.g., Seljak 
2000; Ma & Fry 2000c; Scoccimarro et al. 2001; TJ02; Guzik & Seljak 2002; Sheth et 2001; Takada & Jain 2003a hereafter TJ03a; TJ03b; 
Zehavi et al. 2003). 

Recently, we have extended the halo approach to analytically compute the real-space three-point correlation fvmctions (3PCF) of the 
mass, galaxies and weak lensing fields with reasonable computational expense (TJ03b). The 3PCF is the lowest-order statistical quantity to 
probe non-Gaussianity, generated by non-linear structure formation from primordial Gaussian fluctuations. It also allows one to explore the 
shapes of clustered mass distributions via its configuration dependence, which is not contained in the two-point correlation function (2PCF). 
Hence, measurements of the 3PCF can provide additional cosmological information. Our recent work described above and this paper provide 
the methods and formalism to compute the real space 3PCF. On large scales, the real space 3PCF provides no advantage over the well studied 
bispectrum. On small scales, however, it is easier to measure the 3PCF in real space because the bispectrum requires taking the Fourier 
transform of survey data, which typically involves dealing with the complex survey geometry (e.g. in lensing surveys there are several areas 
that are masked out due to bright stars, and thus the transform is likely noisy). In addition, the different wavenumber modes are highly 
correlated with each other on small scales, and thus one merit of working in Fourier space is lost (in contrast to the CMB, where nonlinearties 
are minimal on scales of interest). Indeed so far most measurements of the 2PCF have also been in real-space for cosmic shear (though Pen, 
Van Waerbeke & Mellier 2002; Brown et al. 2002; Pen et al. 2003 have estimated the shear power spectrum as well). 

In lensing surveys, the 3PCF of the shear is easier to measure than the 3PCF of the convergence field which requires a non-local 
reconstruction from the data (Bemardeau, Van Waerbeke & Mellier 2003). An exciting recent development has been the detection of statistical 
measures based on the shear 3PCF from the Virmos-Descart Survey (Bemardeau, Mellier & Van Waerbeke 2002; Pen et al. 2003). Theoretical 
study of the shear 3PCF has just begun: Schneider & Lombard! (2003, hereafter SL03) and Zaldarriaga & Scoccimarro (2003, hereafter ZS03) 
analyzed how to construct the eight shear 3PCFs from combinations of the +/x shear components at each vertex of a given triangle. Using 
the ray-tracing simulations, in TJ03a we verified that the eight shear 3PCFs display characteristic configuration dependences, and pointed 
out that the characteristics can be used to separate the lensing S-mode from the measured signals that are generally contaminated by the 
S-mode due to systematic errors and other effects. 

The purpose of this paper is to develop an accurate, analytic model to predict the 3PCFs of the shear field, extending the halo model for 
the convergence field (TJ03b). We will carefully examine the accuracy of the model predictions by comparing with ray-tracing simulations. 
In particular, we will focus on whether or not the halo model can reproduce the complex configuration dependences seen in the simulations, 
since the halo model relies on simpUfied assumptions of spherical halos. 

We develop a new application of higher order shear correlations to measure the properties of massive halos. These are relevant for 
upcoming surveys such as the CFHT Legacy Survey^ and the Deep Lens Survey^ and proposed projects such as DML/LSST^, Pan-STARRS* 
and SNAP^. These surveys promise to measure n-point correlation functions of the cosmic shear fields even on sub-arcminutes scales with 
high significance. Therefore, with the halo model, these small-scale signals can be used to probe properties of the halo density profile such as 
its inner slope and concentration. We will work with a generalized universal profile for halos (Navarro, Frenk & White 1996; 1997, hereafter 
NFW). These properties remain uncertain theoretically as well as observationally. To implement this approach, we will develop a method 
to compute the covariance for the 2PCF and 3PCF measurements, extending the method of Schneider et al. (2002b). We will estimate the 
accuracy with which forthcoming surveys can constrain halo profile parameters from combined measurements of the 2PCF and 3PCF. 

The plan of this paper is as follows. In §2, we develop the real-space halo approach for computing the two- and three-point correlation 
functions of the lensing convergence and shear fields. In §3, we summarize the triangle configuration dependences of the shear 3PCFs. In §4, 
we compare the halo model predictions with ray-tracing simulation results. Then, in §5, we show how measurements of the 2PCF and 3PCF 
on sub-arcminute scales are feasible for ongoing and future lensing surveys, taking into account statistical sources of error. We then address 
how these small-scale measurements can be used to constrain halo profile properties. §6 is devoted to a summary and discussion. We will 
consider mainly two cosmological models. One is the ACDM model with Q,mo = 0.3, Q,\o = 0.7, fibo, h = 0.7 and a% = 0.9. The other 
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is the SCDM model with Omo = 1.0, h — 0.5 and as — 0.6. Here Jlmo, ^to and flxo are the present-day density parameters of matter, 
baryons and the cosmological constant, h is the Hubble parameter, and ag is the rms mass fluctuation in a sphere of radius 8/i~^Mpc. 



2 FORMALISM: REAL-SPACE HALO APPROACH TO COSMIC SHEAR STATISTICS 

In this section we develop an analytic method for calculating the 3PCFs of shear fields by extending the real-space halo approach developed 
in TJ03b. 



2.1 Convergence and shear fields 

Weak gravitational lensing can be separated into two effects: magnification (described by the convergence) and shear (e.g., Bartelmann & 
Schneider 2001). 

The lensing convergence field is a scalar quantity and simply expressed as a weighted projection of the density fluctuation field between 

source galaxy and observer: 



= lv'*(0) = r"dxWix)S[x,dA{x)0], (1) 



Jo 



where we have introduced the two-dimensional lensing potential 9, the Laplacian operator defined as V'^ = /dOidOi, x is the 
comoving distance, and xh is the distance to the horizon. Note that x is related to redshift z via the relation dx = dz/H{z), where H{z) 
is the Hubble parameter at epoch z. Following the early work of Blandford et al. (1991), Miralda-Escude (1991) and Kaiser (1992), we used 
two key simplifications to derive the equation above; the flat-sky approximation, which is valid on angular scales of interest, and the Born 
approximation, where the convergence field is computed along the unperturbed path. Using ray-tracing simulations, Jain et al. (2000) showed 
that the Bom approximation is an excellent approximation for the two-point statistics. We will assume that it also holds for the higher-order 
statistics we are interested in. The function W is the lensing projection defined by 

W{x) = lnmoma-\x)dA{x) r"dxs ns{xs) '^^^^'~^\ (2) 



2 7() dA(x 

where ns(xs) is the redshift selection function of source galaxies. Here Hq is the Hubble constant {Ho = lOO/i km s^^Mpc^^) and 
dAix) is the comoving angular diameter distance. In this paper we assume all source galaxies are at a single redshift Zs for simplicity; 
"■s(x) = Sd{x ^ Xs)- Note that dA = X for ^ flat universe. 

A more direct observable of weak lensing is the shearing of images of source galaxies. Since this effect is of order 1% for large-scale 
structures lensing in a CDM model, it is measurable only in a statistical sense. The shear field is described by the two components, 71 and 
72, which correspond to elongation or compressions along the x-axis, or at 45° to it, respectively (given Cartesian coordinates on the sky). 
The shear field is expressed in terms of the lensing potential as 

71 = - *,22), 72 = *, 12, (3) 

where = d^'^ / dOiddj. In Fourier space, these fields are simply related to the convergence field via the relation 

71 (Z) = k{l) cos{2ipi), 72(0 = k{l) sin(2^0, (4) 

where I = l{cos (pi, sin (pi) and, in the following, quantities with tilde symbol denote Fourier components. Equation (4) shows that 7i 
behaves like a spin-2 field: if at a given point one rotates the coordinate system by an angle a in the anti-clockwise direction, the shear fields 
are transformed as 

f[ = cos 2a 71 + sin 2a 72 , 

72 = — sin 2a 71 + cos 2a 72. (5) 

We will often use the vector notation 7 = 71 + 172 . A general two-dimensional spin-2 field can be decomposed into an i5-mode derivable 
from a scalar potential and a pseudo-scalar B-mode (Kamionkowski, Kosowski & Stebbins 1997; Zaldarriaga & Seljak 1997; Hu & White 
1997 for the CMB polarization and Stebbins 1996; Kamionkowski et al. 1998; Crittenden et al. 2002; Schneider et al. 2002a for the cosmic 
shear). Gravitational lensing induces a pure E-mode in the weak lensing regime, while source galaxy clustering, intrinsic alignments and 
observational systematics induce both E and U-modcs in general (Crittenden et al. 2002; Schneider et al. 2002a). The clean separation of 
E/B modes from survey data is of great importance in extracting cosmological information. 



2.2 Halo mass function and halo bias 

To develop the halo model, we begin by describing models of the mass function and the halo bias that are used in this paper. 

Following TJ02 and TJ03b, we employ an analytical fitting formula proposed by Sheth & Tormen (1999), which is more accurate than 
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the original Press-Schechter mass function (Press & Schechter 1974). The number density of halos with mass in the range between M and 

M + dM is given by 

n{M)dM = !^f{y)dv 



^^[l + (az.)--]V^exp(-^)^, (6) 



_ Po , t \-p^ r-T. ( '^^\ '^^ 
_ _^ 

where v is the peak height defined by 

n 2 

(7) 



V 



D(z)a{M) 



Here a{M) is the present-day rms fluctuation in the mass density, smoothed with a top-hat filter of radius Rm = (3M/47rpo)^^^, Sc is the 
threshold overdensity for the spherical collapse model (see Nakamura & Suto 1997 and Henry 2000 for a useful fitting function) and D{z) 
is the growth factor (e.g., Peebles 1980). The numerical coefficients a and p are empirically fitted from iV-body simulations as o = 0.75 
and p — 0.3. Note that the value of a is modified from the a = 0.707 in Sheth & Tormen (1999) to better fit the mass function at cluster 
mass scales in the Hubble volume simulations (R. Sheth; private communication). The coefficient A is set by the normalization condition 
dv}(y) = 1, leading to A = 0.129. This condition reflects the assvunption that all the matter is in halos. Note that the peak hight v is 
specified as a fimction of M at all redshifts once the cosmological model is fixed. 

Recently, various authors (e.g., Jenkins et al. 2001; White 2002; Hu & Kravtsov 2003) have addressed the non-trivial problem of how the 
mass function seen in the simulations depends on the halo identification scheme and the mass estimator. There is no clear boundary between a 
halo and the surrounding large-scale structure, therefore the halo mass does depend on the algorithm used (e.g., the friend-of -friend method or 
the spherical overdensity method). Jenkins et al. (2001) showed that if one employs the halo mass estimator, Miso, enclosed within a sphere 
of radius riso (interior to which the mean density is 180 times the backgrotmd density), the mass function measured from the simulation 
can be well fitted by the universal form of equation (6). This conclusion was also verified by White (2002; also see Hu & Kravtsov 2003). 
Despite these facts, in this paper we employ the virial halo mass to describe the mass function for simplicity, since it is based on the more 
physically-motivated spherical collapse model that can be appUed to any cosmology, irrespective of the halo density profile: 

M = — poA„(0)rei„ (8) 

where po is the present-day mean density of matter, r^ir is the virial radius, and A„(a) is the overdensity of collapse given as a function 
of redshift (e.g., see Nakamura & Suto 1997 and Henry 2000 for a useful fitting formula). One justification of our treatment is the result 
of Figure 5 in White (2002), which showed that the form of equation (6) fits the simulations when the virial mass estimator is employed 
(although the agreement is not as good as the case for Mi8o). The difference in the halo model predictions for lensing statistics due to the two 
definitions is small, as will be shown in Figure A4. Finally, it is worth pointing out the advantage of the mass function of equation (6) over 
that of Jenkins et al. (2001): it is well behaved over the full range of mass and satisfies the normalization condition JJ^ dvfiy) = 1, while 
the mass function of Jenkins et al. (2001) cannot be safely extrapolated outside of the range of their fit and does not satisfy the normalization 
condition. 

Mo & White (1996) developed a useful formula to describe the bias relation between the halo distribution and the tmderljdng mass 
distribution. This idea has been improved by several authors using iV-body simulations (Mo, Jing & White 1997; Jing 1998; Sheth & 
Lemson 1999; Sheth & Tormen 1999); we will use the fitting formula of Sheth & Tormen (1999) for consistency with the mass fimction (6): 

where we have assumed scale-independent bias and neglected the higher order bias fvmctions (62, 63, ■ ■ •). This bias model is used for 
calculations of the 2-halo term in the 2PCF and the 2- and 3-halo terms of the 3PCF. It is not important at the small, non-linear scales where 
the 1-halo term arising from correlations within a single halo provides the dominant contribution. 



2.3 Convergence and shear profiles for an NFW halo 

In this subsection, we derive useful, analytical expressions for the convergence and shear profiles around an NFW halo, which is the most 
essential ingredient for small scales. 

We consider halo density profiles given by the form 

^'^^'"^ " (cr/rvi.)«(l + cr/rvi03-«' 

where ps is the central density parameter and c is the concentration parameter. It is used instead of the scale radius = rvir/c which is the 
transition radius between ph oc r^" and r^'' . For most of the paper, we will use the NFW profile with a = 1. However, since simulations 
with higher spatial resolution have indicated a w —1.5 (Fukushige & Makino 1997; Moore et al. 1998; Jing & Suto 2000; Ghigna et al. 
2000), we will also consider the effect of variations in a for lensing statistics. Given the halo profile, the parameter ps can be eliminated from 
the definition of the virial mass: 

A^pt'- \ f''' (NFW;a = l), 

M= I ^■nr^drphir) = ^^P"^-" { c^"" , . , . (H) 

' ^ ' c3 I !^ 2^1(3 - a, 3 -a, 4- a, -c), (otherwise), 



Jo 
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background galaxy 




Figure 1. Illustration of the shear pattern around an axially-symmetric halo. A background galaxy image is tangentially deformed with respect to the halo 
center, corresponding to a pure E mode in the shear field. 



where / = l/[ln(l + c) — c/(l + c)] and 2^1 denotes the h5rpergeometric function. Note that equation (8) gives the virial radius is given in 
terms of the halo mass M and redshift z. 

To express the halo profile in terms of M and z, we further need the dependence of the concentration parameter c on M and z; however, 
this still remains somewhat uncertain. Following TJ02 and TJ03b we use 

-9 



c{M,z)=co{l + z)-'' 



M 



M,{z = Q) 



(12) 



where M,(z = 0) is the nonlinear mass scale at present defined by 5c{z = 0) = (T(Mt). The redshift dependence {l-\-z)~^ and our fiducial 
choice of (co, /9) = (9, 0.13) are based on the numerical simulation results in Bullock et al. (2001). 

For an NFW profile stated above, we can derive an analytical expression for the the convergence field. As we discussed in TJ03b, the 
halo profile is taken to be truncated at the virial radius in order to maintain mass conservation given by the normalization of the mass function. 
Hence, the convergence field for a halo of mass M can be defined as 



dAiXs) 

where Em is the projected density field defined by: 

Mfc^ 



T.M{e) : 



dri^pM{r\\,dA0) 



2-Rrt 



(13) 



(14) 



The explicit form for F{x) is given by equation (27) in TJ03b. 

For an axially-symmetric mass distribution, the shear amplitude can be expressed in terms of the convergence field as a function of the 
radius from the halo center (e.g., Miralda-Escude 1991): 

'yM{e) = KM{e)-KM{e), (i5) 

where km is the mean surface mass density inside a circle of radius 6: km = {l/nO^) 2n9'd9' km{0'). The equation above reflects the 
non-local nature of the shear field, since it is induced by non-local tidal forces. The shear field does not vanish outside the projected virial 
region because km is non-zero, even if km (6) = 0. From equation (13), the shear profile for an NFW halo of mass M can be analytically 
expressed as 



7m(6I) = A'KGa 
with 



dAiXs) 



(16) 



G{x) = { 



a;2(l + c) 
1 

3(1 + c) 
1 

+ c) 

2f-^ 



(2 - x^)Vc^ - 



l-a;2 
(llc-h 10)Vc2 - 1 



1 + c 

„2\ 



(2 - a;2)Vc2 - a;2 



- 2c 
-6c 
-2c 



2 , 
+ ^ln■ 



c(l+( 



2 - ■ix' 



3.2(1 _ 3.2)3/: 



-arccosh 



x^ + c 
a;(l-^c)' 



2 , 
+ ^In 



x(l + c) 

c + Vc^ -x^ a;2(a;2 - 1)3/2' 



ix-' 



x^ + c 
'a(lTc) 



(a;<l) 
(x = l) 
(1< a: < c) 
{x > c). 



(17) 



where / = l/[ln(l -|- c) — c/(l -|- c)]. Note the asymptotic behavior G{x) — » 1/2 for a; — > 0, while the convergence km — > 00 in this 
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Figure 2. Left panel: Radial profiles of the convergence (solid curve) and sheai' (dashed curve) fields for an NFW halo. The halo mass is A/ = 10 Mq, 
lens redshift zi =0.2 and source redshift Zs = 1. For comparison, the dot-dashed curve shows a singular isothermal sphere {p^ oc r^^) with the same virial 
mass. The two arrows show the scale radius and the virial radius of the NFW profile. Right panel: The angular scales of the scale and virial radii for a lens 
halo are plotted against its redshift. The three curves show the results for halo masses M/Mq = 10^"^, 10^* and 10^^, from bottom to top. Halos relevant 
for lensing statistics have ds ^ 2' for theii' scale radii at redshift z = 0.4. Thus the effect of the inner structures of the halo profile on lensing statistics should 
appear at scales Ji, 2' . 



limit (following from the NFW inner slope ph oc r~^). This feature is in contrast to the case of a power law density profile ph oc r~" 
(1 < n < 3), which leads to km, 7m oc 6^~" and the asymptotic behavior km, Jm oo for 9^0. Outside the virial region the shear 
amplitude behaves like the field around a point mass, "fm oc 1/9^, due to the sharp cutoff of the mass distribution. Taking c — > cxa in the 
equation above leads to the expression for 7m in Bartelmann (1996; also see Wright & Brainerd 2000), which is derived from a line-of-sight 
projection of the NFW profile under the assumption that the profile is valid for an infinite range beyond the virial region. As shown in Figure 
A4, if one adopts the expression for jm in Bartelmann (1996) the amplitude of shear correlations is significantly overestimated. Therefore, 
the halo profile boundary should be carefully considered to obtain accurate results as well as a consistent formulation. 

The sketch in Figure 1 illustrates how a background galaxy image is tangentially deformed around a foreground lens with an axially- 
symmetric profile. In the weak lensing regime, the shear pattern is described by a pure _E-mode. In Cartesian coordinates with the center 
taken as the halo center, the two shear components can be expressed as 

7m,i(0) = -7Af(6') cos2(^0, lM,2{d) = -7Af(6') sin2(/3g, (18) 

where tpg is the angle between the first x-axis and the line connecting the halo center and the source galaxy position (see Figure 1). 

The left plot in Figure 2 shows the radial profiles of the convergence (solid curve) and shear (dashed curve) for an NFW halo of mass 
10^'^h~^ Mq, computed using equations (13) and (16). We use the ACDM model, lens redshift zi — 0.2 and source redshift Zs = 1. The 
two arrows denote the angular scales 9b and 5vir corresponding to the scale radius and the virial radius, respectively. In contrast to the 
convergence, the shear is almost constant in the inner region and does not vanish for 9 > 9vir, rather it follows the point mass profile cc 9~^. 
It should be noted that the shear profile appears to be not a smoothly varying function at ^ = 9vir because of the artificial sharp halo boundary. 
For comparison, the dot-dashed curve shows the shear amplitude for a singular isothermal sphere (SIS) with the profile ph oc r^^ (for this 
case the convergence and shear fields coincide). The comparison manifests the characteristic feature for the NFW shear field. For example, 
the plateau profile, 7m ~ constant, at 9 ^ 9s could be a direct test of the inner slope ph oc r~^. 

The right panel in Figure 2 explicitly shows the angular scales of the scale and virial radii for lens halos as a function of redshift. We 
consider the mass scales, M/Mq = 10^^, 10^^ and lO'^^. Halos with M > W^^Mq, which dominate the contribution to lensing statistics 
on small, non-linear scales, have 9s ^ 2' for their scale radii at redshift z — 0.4. Thus, the figure implies that the the inner region likely 
affects the lensing statistics on angular scales ^ 2', as we will show explicitly in Figure 13. 



2.4 Real-space halo approach 

In the halo model, the 2PCF of the convergence field, ^k{9), can be expressed as the sum of correlations within a single halo (1-halo term) 
and between two halos (2-halo term). In TJ03b, we obtained the real-space relation for the 1-halo term: 

(e) = M0i)^(02))^'^ = J'''dx-^JdMn{M;x)J^^^dsJ dip s km{s;x)^m{\s + e\;x), (19) 
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Figure 3. Projected number density of halos more massive than a given mass scale M, N{> M), between z = and 1 for the ACDM model. 



where cfV/dxd^l = d\{x) for a flat universe, and we have used polar coordinate s = s(cos (p, sin tp) to write d^s = sdsdip. {■ ■ ■) 
denotes the ensemble average. From statistical symmetry, we can set the separation vector 6 to be along the first axis, so that \s + d\ = 
(s^ + 2se cos (p)^/^ To obtain ^l'' for a given cosmological model, we perform a 4-dimensional numerical integration. Equation (19) 
shows that ^i'' is given by the sum of lensing contributions due to halos along the line of sight weighted with the halo number density on the 
light cone. Hence, this form correctly accounts for multiple lensing due to halos at different redshifts. 

Figure 3 plots the angular number density of halos more massive than a given mass scale M between the source redshift 2s = 1 and 
the present: iV(> M) = J^'dz d^V/{dzdn) J^dM' n(M'; z). One can see iV(> M) < 10"^ arcmin"^ for halos with M > lO^^M©, 
which provide dominant contribution to the lensing statistics on small angular scales (see Figure 14). This clarifies that the 1-halo contribution 
is primarily due to a single lens halo, and that there is only a small probability for multiple lensing due to such massive halos at different 
redshifts. However, one should keep mind the importance of multiple lensing due to such cluster-scale halos and less massive halos, as shown 
by White, Van Waerbeke & Mackey (2002) and Padmanabhan, Seljak & Pen (2003). 

The form of equation (19) allows us to straightforwardly extend it to compute the 2PCF of shear fields by replacing the convergence 
profile. Km, with the shear profile, 7m, for a given halo. The 1-halo term in Cartesian coordinates is 

(7^(^1)74^2))"'=/ dxd\{x) IdMniM-x) ds d^ s -iM{s;xhM{\s + e\■x)e^.{s)e.{s + 0), 
Jo J Jo Jo 

where is the phase factor of the shear field as given by equation (5); e/j(s) = —(cos 2^5, sin 2^5) and ep(s + 0) = 
— (cos 2(^g_l_g, sin 2(/3g^g) with {cos ipg^Q, sirup ^^q) = —\s + 0\^^{scosp + 9,ssinp). In contrast to equation (19), the equation 
above employs an infinite integration range for d^s in order to account for the non-local property of the shear fields. In practice, setting the 
upper bound of Jds to be three times the projected virial radius gives the same result, to within a few percent. We investigate the accuracy of 
equation (20) as well as its self-consistency with other methods in § A. 

Since the two shear components are not invariant under coordinate rotation, the shear 2PCFs generally depend on the relative orientation 
of the two points (e.g.. Kaiser 1992). This issue has been well studied in the literature (Kamionkowski et al. 1998; Crittenden et al. 2002; 
Schneider et al. 2002a). One way to avoid the coordinate dependence is to use the -|-/x decompositions, where the + component is defined 
as the shear field in parallel or perpendicular direction relative to the Une connecting the two points taken, while the x component is the 45° 
rotated shear field. We can thus define the rotationally invariant 2PCFs of the shear field: 



(20) 



MO) = M0ih+(02)), 

C7.x(e) = (7x(0l)7x(02)). 



(21) 



It should be noted that (7+(0i)7x (62)) and (7x (0i)7+(02)) vanish because of invariance under parity transformation. We will also 
consider the following 2PCF: 



e4e) = (7(0i)-7*(02)). 



(22) 



Note that ^^(0) = ^7(^)- The 1-halo term contributions to these shear 2PCFs can be calculated by replacing the phase factors €i^e„ in 
equation (20) with 



cos 2^5 cos 2(pg^g -I- sin2(^sin2(^g_l_g, for ^-y, 
ef,{s)e^{s + 0) ^ { cos2(^cos2(/3g_|_^, for 5^,+ , 

sm2ipsm2ipg^g, for f^,x, 



(23) 
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mirror transformation 




Figure 4. Upper. A sketch of the triangle configuration, parameterized by r, q and ip, used to describe the 3PCF. The solid and dashed lines at each vertex 
show the positive directions of the + and X components of the shear field, defined with respect to the triangle center denoted by o. Lower. The mirror 
transformation of the triangle with respect to the side vector r for $x++ is shown as one example. The 7x component at each vertex changes sign under the 
mirror transformation (at vertices 1 and 1' in this case). 

An alternative approach to the lensing 2PCF, which is conventionally used in the literature, is based on a model of the non-linear 3D 
power spectrum of the mass, P{k). Combining P{k) with Limber's approximation (Limber 1954; Kaiser 1992) allows us to compute the 
shear 2PCFs: 

(7m7.>W = j^'dx W\x,Xs)d-/{x) P{k = File). (24) 

Setting the window function to F{x) = Jo{x), [Jo{x) + J4{x)]/2 and [J()(x) — J4(x)]/2 yields and (,-,,x, respectively. So far, 

there are two well studied models for the nonlinear P{k). One is the fitting formula calibrated from iV-body simulations (e.g., Jain, Mo 
& White 1995, Peacock & Dodds 1996, and Smith et al. 2002). This method has been extensively used for the interpretation of cosmic 
shear measurements in terms of cosmological parameters (e.g.. Van Waerbeke 2001b). The other is the recently developed Fourier-space halo 
approach (Seljak 2000; Ma & Fry 2000b,c; Scoccimarro et al. 2001). In the halo model, P{k) is similarly expressed as sum of the 1- and 
2-halo terms; P{k) = Pih{k) + P2h{k) (see §2.2 in TJOb for details). Inserting P\h{k) into equation (24) leads to the 1-halo term of the 
shear 2PCFs in the Fourier-space halo model. Note that the condition ^„ = holds, since they both have the same window function Jo in 
equation (24). 

We turn to the 3PCF of lensing fields, which is main focus of this paper. The 1-halo term of the 3PCF of the convergence field is given 
by equation (52) in TJ03b: 

C^''(r,g,V) = (k(6'i)k(6»2)k(6»3))"' = / dx d\{x) dM n{M;x) ds d^p s km{s)km{\s + r\)KM{\s + q\), (25) 

Jo J Jo Jo 

where we have set |s + r| = (s^ + + 2rscos'fiY^^ and |s + g| = [s^ + <f + 2sqcos{(p — tjj)]^^^. From statistical symmetry, the 
convergence 3PCF can be expressed as a function of three parameters r, q and characterizing the triangle configuration, as shown in Figure 
4. Note that we often omit x in the argument of km or jm for simplicity. This equation means that we can compute (^^.^ by a 4-dimensional 
integration for any triangle configuration, which is the same level of computation as the 2PCF. This holds for higher-order moments as well 
(TJ03b), which is a great advantage of the real-space halo model. For comparison, the Fourier-space halo model requires a 6-dimensional 
integration to get Ci^. In TJ03b, we have also developed approximations for computing the 2- and 3-halo terms of Ck. The 2-halo term is 
relevant for a range of transition scales between the non-linear and quasi-nonlinear regimes, while the 3-halo term dominates in the large 
scale, quasi-linear regime where perturbation theory (PT) is valid. 

We can extend equation (25) to the shear 3PCFs. Since the shear field has two components at each vertex of a given triangle configuration, 
we can formally construct eight (2^ = 8) components of the shear 3PCFs. In Cartesian coordinates the functions are 
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C-y,abc{Ol,02, Os) = {l a{e l)')i{e 2)1 dO 3)) , (26) 

where a,b,c — 1, 2. Note that the subscripts a, b, c correspond to the shear components at vertices 9i, O2 and 63, respectively, in our 
convention (see Figure 4). As for the shear 2PCFs, the 3PCFs defined above are not invariant under coordinate rotation. Therefore, in contrast 
to the 3PCF of the convergence or any scalar quantity, they are not specified by three parameters. Instead four parameters are needed, e.g., 
r, q, ij) and the orientation angle of the side vector r relative to the first coordinate axis. The real-space halo model yields the following 
expressions for the 1-halo terms of the shear 3PCFs: 

C^,ahckr,q) = \ dxdi(x) dM n{M;x) ds dip s 'yM{s)^M{\s + r\)^M{\s + q\)ea{s)eb{s + r)ec{s + q), (27) 
Jo J Jo Jo 

with 

ea(s) = — (cos 2<^, sin2(^), 
eb{s + r) = — (cos2</9s+r,sin2((9s+r), 

ec(s + q) = -{cos2ips+q,sin2(fis+q) (28) 

where {cos (ps +r, sirups +r) — (r + s cos tp, s sin (p)/\s + r\, (cos (ps+q, sin ips+q) = (gcos?/' + scosip,qsinip + ssintp)/\s + q\, 
\s + r \ and |s + q| are given below equation (25), and we again take the halo center as the coordinate center, using statistical symmetry. 

To avoid the coordinate dependence for the shear 3PCFs, we again use the +/x decomposition of the shear fields, as stressed in SL03. 
For the three-point case, however, there is no unique choice of the reference direction to define the -|-/x components. Following ZS03 and 
TJ03a, we take the 'center of mass', o, of the triangle, defined by 

o = -''^Oi = s + -{r + qcosip,qsinip). (29) 
The projection operators, which transform 71 and 72 at each vertex into the +/x components, with respect to o, are 

p+{e,) = -{el - el, 2ene,2)/el Px {0,) = -(-2^afe, el - el)/el (30) 

where Oi = Oi — o (i = 1,2, 3). The geometry of the triangle we consider is illustrated in Figure 4, where the solid and dashed lines at each 
vertex denote positive directions of the -|- and x components. Using these projections, we can thus define the shear 3PCFs from combinations 
of the -I-/ X components so that the resulting 3PCFs are invariant under rotations of the triangle with respect to o: 

C7..Mr('',?,V') = P;(0l)P.^02)P=(03)(7a(0l)76(02)7c(03)>, (31) 

where /i, z/, r = + or x . Inserting equation (27) into the r.h.s above yields the 1-halo term for C,~f^v^T- The projection operators and (• • •) 
commute, since the projection operators do not depend on the integration variable s, as seen from equation (30). The shear 3PCFs defined 
in this way are functions of the three parameters r, q and ■)/'• Although we adopt the center of mass throughout this paper, SL03 proved that 
the eight shear 3PCFs with respect to any center can be expressed as linear combinations of the eight 3PCFs above. To complete the halo 
model predictions, it is necessary to develop the perturbation theory predictions for the shear 3PCFs, which are relevant on large scales. For 
the convergence we have obtained these, but the complex spin-2 properties of the shear fields make it non-trivial to obtain the PT prediction. 
We will consider only the 1-halo term for predictions of the shear 3PCFs. Nevertheless, the 1-halo term agrees with the simulation results on 
scales of interest, as shown below. 



3 TRIANGLE CONFIGURATION DEPENDENCES OF THE SHEAR 3PCF 

From the observation that there are eight shear 3PCFs, SL03 and ZS03 addressed the questions: how many functions are non-zero? How does 
each function carry information about the E/B-modesl In TJ03a, using ray-tracing simulations, we qualitatively verified the conclusions of 
SL03 and ZS03, which were based on analytical studies. We found that: (1) The eight 3PCFs are generally non-zero. (2) For a pure i?-mode, 
two or four components vanish for isosceles or equilateral triangles, respectively. (3) We also studied the triangle configuration dependence 
of the shear 3PCF given by equation (36), and pointed out that they offer a promising way to disentangle the E/B-modes for measured shear 
3PCFs. Here we summarize these characteristics of the shear 3PCFs. 

We restrict ourselves to a pure E field, expected in the weak lensing regime. We consider a mirror transformation as shown in Figure 4, 
which illustrates the transformation with respect to the side vector r for Cx++- It corresponds to t/; — > 27r — -tp in our parameterization. From 
statistical homogeneity and symmetry, the amplitude of the shear 3PCF depends only on the distances between the center and each vertex. 
Hence, the absolute amplitudes of Cx++ for the two triangles shown should be same. But the sign of 7x at the vertex 1' changes under this 
mirror transformation. According to this property, we can divide the eight 3PCFs into two groups: 

Parity-even functions: CnuT{r,q,ip) = CiJ.vr{r,q,27r - ip), for (/x, i^,r) = (-h, -h, -h), (-h, x, x), (x, -h, x), (x, x, -h), 
Parity-odd functions: C^^r{r,q,ip) = -C^>„T{r,q,2^^ - ip), for {ij.,i',t) = {x,x,x),{x, +,+),{+, x, +),{+, +,x). (32) 

Note that the 3PCF of a scalar quantity transforms as (^(r, q, tp) = (^(r, q, 2n — ■ib). 

Next we consider special triangle configurations: the first is an isosceles triangle with r = g. In this case, the 7x components at 
vertices 1 and 1' in Figure 4 are statistically identical (viewed from the center of the triangle, they should have equal contributions when 

© 0000 RAS, MNRAS 000, 1-31 



10 M. Takada & B. Jain 



averaged over the matter distribution)^. We thus have additional symmetries for the two 3PCFs: ^x++('', q, tp) = Cx++{r, q, 27r — ip) and 
Cx X X (2^1 ,xi,ip) — (^xxx{xi,xi,2tt — tp). These relations and equation (32) yield 

Isosceles Triangles: Cx++ = Cx x x =0. (33) 

Note that the other two parity-odd functions, C+x+ and C++x , do not vanish in general, since the component 7x is at a vertex bounded by 
unequal sides. For equilateral triangles, however, all four parity-odd functions vanish: 

Equilateral Triangles : Cx x x = Cx++ = C+x+ = C++x = 0. (34) 

The 3PCFs discussed above do not vanishing for a B-mode shear field, as explained below. Hence these 3PCFs provide a direct, simple test 
of the J3-mode contribution to the measured signal. 

Let us consider again generic triangle configurations, but now for a pure B field. As discussed by TJ03a, the analog of equation (32) for 
a B-mode spin-2 field is: 

C^^r{r,q,tl}) = -CMi'r(r,g,27r-V), for {^i,u,T) = (+, -h, -h), (+, x , x ), (x , +, x ), (x , x , +), 

CM-'r(r,g,?/)) = Cf,^r{r,q,2TT-i>), for {^i,v,T) = (x , x , x), (x ,+,+),{+, x ,+),(+, +, x). (35) 

Thus the symmetric and anti-symmetric functions are reversed compared to the i5-mode. This follows from the fact that given a pure E field, 
we can generate a pure B field by rotating the E field at each point by 45 degrees (Kaiser 1992). We can then consider the eight 3PCFs of 
a pure B mode similarly as done for the shear 3PCFs. Since this procedure transforms the original S-mode components 7+ and 7^ at each 
point into 7^ and 7^ in the transformed B field, respectively, we get C+++ Cx x x and so on. For a general spin-2 field that contains both 
E/B modes, the configuration dependences of equations (32) and (35) do not hold. The eight 3PCFs therefore have to be measured over the 
full range i/> = [0, 2tt], unlike the 3PCF of a scalar quantity. 

From the symmetry properties discussed above, we propose simple estimators for the E/ B-mode contributions to measured shear 
3PCFs: 

Estimator of i;-mode : C^r{r,q,'4') = ^ [C-f,i^,^T{r,q,ip) ± C-(,P"'r{r,q,2Tt - ijj)], 

Estimator of S-mode : C^^^ 9, V') = ^ [Ct.m-^t (r, g, ■0) T Ct.m-'t- {r,q,2Tr - ip)], (36) 

where the upper and lower signs in ± or =p are meant for [ii,v,t) = (-|-, -|-, -I-), (+, x , x), (x , +, x ), (x , x , -I-) and 
(x, X, x), (x, -|-, -|-), X, -I-), x), respectively. Note that this argument holds if the E and B modes are statistically uncorre- 

cted. These estimators allow one to separate the lensing S-mode contribution from the measured 3PCFs that are in general contaminated by 
the _B-modes contribution of intrinsic alignments, source galaxy clustering and observational systematics. In comparison, for the case of the 
shear 2PCFs, a non-local integration is required to discriminate the lensing iJ-mode (Schneider et al. 1998; Crittenden et al. 2002; Schneider 
et al. 2002a). 

After the subirussion of this paper, Schneider (2003) pointed out that if a i? field is parity invariant in a statistical sense, any correlation 
function that contains an odd number of B-mode shear components vanishes. This is likely to be true for a cosmological B field such as 
intrinsic alignments if wc have a sufficient survey area. Therefore, the argument wc have made above is valid only for a parity non-invariant 
B field. This can be seen in that the 45 degree rotation procedure described above to generate a B field from simulations would produce shear 
fields around clusters with a clockwise curl direction, which violates statistical parity invariance. Nevertheless, we believe our discussion 
is useful, because generic observational systematics are likely lead to the violation of parity invariance. Hence, these arguments strength 
practical usefulness of the shear 3-point functions to disentangle E/B modes from the measurement. 



4 COMPARISON WITH RAY-TRACING SIMULATIONS 

In this section, we address the accuracy of the halo model for predicting lensing statistics, in particular the 3PCFs, by comparing model 
predictions with ray-tracing simulation results. The validity of the real-space halo model for shear correlations, developed in §2, is carefully 
investigated in §A. 

4.1 Cosmological models 

In this paper, we mainly consider two CDM models whose cosmological parameters are chosen to facilitate comparison with the ray-tracing 
simulations used below (Jain et al. 2000; Menard et al. 2003; Hamana et al. 2003). One is the SCDM model {Q,mO = 1, h — 0.5 and 
0-8 = 0.6), and the other is the ACDM model (Qmn = 0.3, Oao = 0.7, Qto = 0.04, h = 0.7 and as = 0.9). For the ACDM model, we 
need to care about the baryon contribution to the input primordial power spectrum. Although we assume a scale-invariant power spectrum 
for the primordial fluctuations, we employ different CDM transfer functions for the SCDM and ACDM models. For the SCDM model, we 
use the transfer function in Bond & Efstathiou (1984) with the shape parameter F = Qmoh = 0.5. On the other hand, for the ACDM model 
we employ the BBKS transfer function (Bardeen et al. 1986) with the shape parameter in Sugiyama (1995), since the shape parameter 
approximately describes the baryon contribution. 

^ This argument is true only for a pure E field, since it relies on the invariance of the B-mode under parity transformation. 
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4.2 Ray-tracing simulations 

We use ray-tracing simulations of tlie lensing convergence and shear fields. We will mainly consider two cosmological model simulations, 
the ACDM model (kindly made available to us by T. Hamana; see Menard et al. 2003; Hamana et al. 2003 for details) and the SCDM model 
(Jain et al. 2000), as described in §4.1. The A'^-body simulations on which these are based were carried out by the Virgo Consortium ^ (also 
see Yoshida, Sheth & Diaferio 2001), and were run using the particle-particle/particle-mesh (P^M) code with a force softening length of 
'soft ~ 30/i~^kpc. The linear power spectra of the initial conditions of the simulations were set up using the transfer function from CMBFast 
(Seljak & Zaldarriaga 1996) for the ACDM model and the fitting function in Bond & Efstathiou (1984) for the SCDM model, respectively. 
For the halo model predictions for ACDM , we employ the BBKS plus Sugiyama transfer function as described above. We verify in Figure 
A3 that the transfer function is sufficiently accurate (see also Eisenstein & Hu 1999). 

The ACDM simulation employs 512^ CDM particles in a cubic box 479/i~^Mpc on a side, while the SCDM simulation uses 256^ 
particles in a 84.5/i~^Mpc side-length box. The particle mass is rripart = 6.8 x IQ^^hT^ Mq and 1.0 x lG^^h~^ Mq, respectively. The 
resolution scale, which is not affected by the discreteness of the A'^-body simulations, is estimated roughly as being ten times the grid size, 
leading to A 94/i~^kpc and SS/i^'^kpc for the ACDM and SCDM models, respectively. Therefore, the SCDM model simulation has better 
spatial resolution than the ACDM model. 

The multiple-lens plane algorithm to simulate the lensing maps from the A/^-body simulations is detailed in Jain et al. (2000) and Hamana 
& Mellier (2001). Throughout this paper, we use a single source redshift of Zs = 1. The simulated areas of the lensing maps are Q,s = 11.7 
and 7.69 degree'^ for the ACDM and SCDM models, respectively. The map is given on 1024^ grids with grid spacing ^grid = 0'.2 and 0'.16 
for the two models. To analyze the correlation functions we should be careful about the possibility that the finite angular resolution could 
affect computations of the n-point correlation fimctions from the simulated maps on small scales. It is not easy to infer the effective angular 
resolution from the spatial resolution of the A'^-body simulations due to the broad lensing projection kernel. Further, the projected density 
field was smoothed to suppress the discreteness effect of the A^-body simulations. The smoothing scale is likely to determine the angular 
resolution rather than the spatial resolution of the A^-body simulation. Note that the ACDM simulation we employ below is the one labelled 
with small-scale smoothing in Menard et al. (2003). These resolution issues were carefully discussed in Menard et al. (2003) and Jain et al. 
(2000), as can be seen from Figure 3 in Menard et al (2003) and Figure 2 in Jain et al. (2000). The former shows the projected angular scale 
of the smoothing as a function of redshift, while the latter shows the angular scale of the force-softening length of the A^-body simulation. 
These scales are ^rcs ~ 0'.3 and 0'.2 at z = 0.4 for the ACDM and SCDM models, respectively, where z — 0.4 is approximately the peak 
redshift of the lensing efficiency for source redshift Zs = 1. Therefore, angular scales that are not affected by finite resolution are likely to 
be ^1', although the resolution of the ACDM model simulation might be slightly worse than that of the SCDM model, as stated above. We 
will keep in mind these resolution issues in the following analysis. 

To compute the sample variance of the lensing statistics from the simulations, we use 36 and 9 realizations of the simulated lensing maps 
for the ACDM and SCDM models, respectively. The realizations are generated by randomly rotating and translating the N-body simulation 
boxes (using the periodic boundary conditions of the A^-body boxes) when the ray-tracing simulations are performed. However, they were 
built from one realization of the A/^-body simulation, which is a sequence of the redshift-space evolution of large-scale structure. Therefore 
the different realizations of the simulated lensing maps are not fully independent of each other. In particular, this could matter when we 
compute the covariance for the n-point correlations in the different bins from the simulations. This issue is still an open question, to be 
addressed in the future using a sufficient number of truly independent ray-tracing simulations. 

4.3 Algoritlim for computing the 3PCF from simulated maps 

To compute the 3PCFs of the lensing fields from the ray-tracing simulations, we implement the method described in §3.2.2 in Barriga & 
Gaztaiiaga (2002). The 3PCF is given as a function of three parameters (r, q and tp) specifying the triangle configuration (see Figure 4). 
The question is how we can efficiently find triplets from the simulated lensing map with Agrid grid points, subject to the constraint that the 
triplet forms a given triangle configuration within the bin widths. We first select vertex 1 on the grid, and then search vertex 2 in the upper 
half plane in an annulus of radius r with given bin width, centered on vertex 1. For given pair 1-2, we look for vertex 3 in a semi-annulus 
of radius q in the upper plane above the line connecting vertices 1 and 2 (in the anticlockwise direction from vertex 1 as our convention), 
imposing the condition that the three vertices form the required triangle configuration within the bin widths. This results in all triangles being 
counted once, if we impose the conditions 9i2 < 623 < 631 . We can use the same list of neighbors to find vertices 2 and 3 for each vertex 1. 
This process is illustrated in Figure 3 in Barriga & Gaztafiaga (2002). Further, to compute the 3PCFs of the shear fields, we need to compute 
the -h/x components of the shear fields at each vertex of the triangle. The projection operators to compute these components are given as a 
function of the list of neighbors of vertices 2 and 3, independent of the position of vertex 1, as can be seen from the definition of equation 
(30). In summary, the 3PCF computation from the simulation map requires roughly O(Agrid) operations on sufficiently small scales. This is 
significantly faster than a naive, direct implementation which requires 0{N^^i^) operations. 

4.4 The two-point correlation function of the shear fields 

In Figure 5 we present a comparison of the halo model predictions for the shear 2PCFs with those measured from the ray-tracing simulations 

for the ACDM (upper panel) and SCDM (lower panel) models. The left panel is the comparison for The solid curve is the halo model 
prediction, while the square symbol denotes the simulation result. The dashed curve denotes the result computed from the Smith02 formula. 

see http : / /star-www. dur . ac .uk/ "f razerp/virgo/virgo . html for the details 
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Figure 5. Comparison of model predictions for the shear 2PCFs with measurements from the ray-tracing simulations. The left panel shows the shear 2PCF, 
= (t ■ T*), for the ACDM (upper panel) and SCDM (lower panel) models. The solid and dashed curves are the predictions from the halo model and the 
Smith02 formula, respectively, while the square symbols with eiTor bars denote the simulation result. Note that errors in different bins are highly correlated 
with each other Similarly, the right panel shows the comparison for §7,+ ™d x ■ 



The theoretical predictions agree very well with the simulation results for the two CDM models. The smallest scale data from the simulations 
corresponds to three times the grid size, and the error bar in each bin is the sample variance for area Q, = 11.7 and 7.7 degree'^ for the 
ACDM and SCDM model simulations, respectively, as estimated from the scatter among their 36 and 9 realizations (see §4.2). We note that 
the errors in different bins are highly correlated with each other. Even if one combines two neighboring bins, the error amplitude remains 
almost unchanged. We have confirmed that this is true for all the statistical quantities we consider below, over angular scales of interest. 

Similarly, the right panel in Figure 5 shows the comparison for ^+ and ■ The theoretical predictions for are again in agreement 
with the simulation results. On the other hand, for (,-y,x, the halo model prediction lies slightly below the simulation results at 6* ^ 5', while 
the Smith02 prediction agrees better. If we employ the halo boundary of rigo as discussed in Figure A4, the model predictions agree better 
with the simulations. The relation between the amplitudes of ^^,4. and ^^^x is physically determined by the fc-slope of the underlying mass 
power spectrum, as can be seen from equation (24). Within the framework of the halo model, the slope is determined by the combined effects 
of the halo profile, the mass function and the slope of the primordial power spectrum (Seljak 2000; Ma & Fry 2000a,b,c; Scoccimarro et al. 
2001; TJ02; TJ03b). Thus, separate measurements of ^+ and ^x could constrain these physical ingredients. 



4.5 The three-point correlation function of the convergence field 

In the following, we address the accuracy of halo model predictions for the 3PCFs of lensing fields by comparison with ray-tracing sim- 
ulations. Until our recent work (TJ03a,b), there had been no analytical model for the lensing 3PCFs on small angular scales, except for 
investigations of the skewness and bispectrum of the convergence field based on extended perturbation theory (Scoccimarro & Frieman 
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Figure 6. The 3PCF of the convergence field for equilateral triangles against side length r in arcminutes. Top panel: The thick solid curve shows the halo 
model prediction for the reduced 3PCF, Qk., defined by equation (37). The thin solid curves denote the 1-, 2- and 3-halo term contributions separately. The 
dashed curve is the perturbation theory (PT) prediction. Middle panel: The comparison with simulations for Qk., as in Figure 5. The simulation result does not 
display the bump in Qk. seen in the halo model prediction. This is likely to be due to halo boundary effects in the standard implementation of the halo model 
(see text for details). The dot-dashed curve denotes the result if we replace the 2- and 3-halo terms with the PT prediction. This modified halo model agrees 
well with the simulation result. Bottom panel: The comparison for the 3PCF itself, as in the middle panel. 



1999; Hui 1999; Van Waerbeke et al. 2001a) or the Fourier-space halo model (Cooray & Hu 2001a; TJ02). Hence we present a detailed 
analysis of the accuracy of halo model predictions for the lensing 3PCFs. 

First, we consider the 3PCF of the convergence field. As in the literature, we define the reduced 3PCF as 

^ , , N C,K{r,q,ip) 

^KKr,q,V) - ^^i^r)UD+Ur)U\r - q\) + Uq)U\r - q\) ' 

where we have used the parameters r, q and i/' to describe the triangle configuration and \r — q\ — (r^ + q^ — 2rq cosi/;)^^^ (see Figure 4). 
The reduced 3PCF is sensitive to QmO, but insensitive to the power spectrum normalization ag, scaling roughly as Qk oc f2~o (Bernardeau 
et al. 1997; Jain & Seljak 1997; TJ02). Therefore, measuring the 3PCF is expected to break degeneracies in the determination of as and Slmo 
from measurements of the shear 2PCFs. 

Figure 6 shows the convergence 3PCF for equilateral triangles against side length for the ACDM model. The thick solid curve in the 
top panel shows the halo model prediction for Qk- A bump feature is evident over the range 1 ^ r ^ 5'. As discussed below, this feature is 
unlikely to be real. The three thin curves show the 1-, 2- and 3-halo terms separately. These contributions to the convergence 3PCF are present 
in the numerator of Qk, but the 2PCF in the denominator includes the full contribution (1- plus 2-halo terms). The 1-halo term provides the 
dominant contribution to Qk at small scales, r ^ 3'. The 2-halo term is relevant over the transition scales between the non-linear and linear 
regimes. The bump feature in Qk is mainly due to the 2-halo term contribution. The 3-halo term eventually dominates on larger scales, and 
gives the perturbation theory (PT) result for Qk at r ^ 10'. The 2- and 3-halo terms appear to be relevant over a wide range of angular scales 
compared to the 3D mass 3PCF, Q, shown in Figure 7 in TJ03b, since the lensing projection causes various length scales at different redshifts 
to contribute to the lensing statistics. Even at the non-linear scale of r = 1', the 2- and 3-halo terms make 11% and 4% contributions to Qk- 

The middle panel shows a comparison of the halo model prediction with the simulation result for Qk, as in Figure 5. The bump feature 
in the halo model prediction cannot be seen in the simulation result and, as a result, the halo model overestimates the simulation result at 
more than la. We have confirmed that this is also true for the SCDM model (also see Figures 8 and 10 in TJ03b for the discrepancy between 
the halo model prediction and the simulation result for the 3D mass 3PCF). There are two effects to be considered in finding the origin of 
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Figure 7. Plots for the SCDM model, as the middle panel in Figure 6. The square symbols show the simulation result for for the SCDM model. The 
triangles denote the simulation results for the ACDM model, with the error bars scaled to the same area. The comparison shows the sensitivity of Qk, to the 
cosmological model, in particular to fi^o ■ 

the bump feature. First, the standard halo model does not take into account halo exclusion effects: the 2- and 3-halo terms should require 
that different halos be separated by at least the sum of their virial radii. Second, so far we have used a sharp cutoff for the halo profile, 
which could lead to inaccuracy in the prediction as discussed for Figure A4 In fact. Figure 8 in TJ03b clarified that modifications aimed at 
resolving these two effects do suppress the bump feature in the 3D mass Q (see also Somerville et al. 2001, Bullock et al. 2002 and Zehavi 
et al. 2003 for discussions on the halo exclusion effect). However, resolving these problems requires a more careful and systematic study, in 
combination with A^-body simulations. This is beyond the scope of this paper. Therefore, we instead employ a simple prescription to avoid 
the overestimation of ~ we replace the contribution from the 2- plus 3-halo terms in the 3PCF with the PT prediction. This treatment 
preserves two merits of the halo model: the quasi-linear 3PCF at large scales is reproduced by the PT result alone, and the non-linear 3PCF 
on small scales is described primarily by the 1-halo term. The dot-dashed curve shows the modified halo model prediction for Qk - It displays 
excellent agreement with the simulation results over the scales we have considered. An alternative method is to simply ignore the 2-halo term 
contribution to the 3PCF, leading to similar agreement. Hereafter, we will use the 1-halo term plus the PT result as the halo model prediction 
for the convergence 3PCF. 

The halo model and the simulations show a flattening of at small scales, r 3' (the halo model predicts a slight decrease of Qk 
with decreasing r). The same feature is found for the SCDM model, as shown in Figure 7. This is consistent with the hierarchical ansatz for 
relating higher order correlations with the 2PCF (e.g. Peebles 1980). For the halo model this behavior results mainly from the NFW profile 
and the halo concentration used (Bullock et al. 2001). It will be of great interest to address these issues in detail using ray-tracing simulations 
with higher angular resolution that probes sub-arcminute scales. 

The bottom panel compares the halo model prediction for the 3PCF (1-halo term plus PT) with the simulation result. The agreement 
is excellent, while it is clear that PT substantially underestimates the 3PCF at scales 9^5'. However, in contrast to the middle panel, the 
standard halo model prediction (lh-l-2h-l-3h terms) denoted by the solid curve displays equally good agreement. This explains why the 
parameter is a more sensitive indicator of gravitational clustering (see Bemardeau et al. 2002b for an extensive review). For example, the 
amplitude and configuration dependence of the Q parameter have been widely used in the literature in connection with questions of stable 
clustering and the hierarchical ansatz (e.g., Peebles 1980; also see Jain 1997; Ma & Fry 2000b; TJ03b). In fact, from a comparison between 
the middle and bottom panels, one can see that the Qk parameter displays a pronounced transition between the non-linear and quasi-linear 
regimes - Qk ~ 45 and ^ 30 at S ^ 3' and ^ 5', respectively. This feature originates from the transitions in the 2PCF and 3PCF of the 3D 
mass distribution predicted for CDM structure formation (e.g.. Figures 1 and 1 1 in TJ03b). 

Figure 7 shows a comparison of the halo model prediction for Qk (dot-dashed curve) with simulation results (square symbol) for the 
SCDM model, as in the middle panel of the previous figure. The halo model prediction matches the simulation results. For comparison, the 
triangle symbols denote the simulation results for the ACDM model, which shows the strong sensitivity of Qk to the cosmological model, 
especially to fimo. 

The accuracy of the halo model is further explored in Figure 8, which compares the prediction for the reduced convergence 3PCF, Qk , 
with simulation results against triangle configurations for the ACDM model. Since the halo model employs the spherically symmetric NFW 
profile, it is interesting to examine whether or not the halo model can properly describe the configuration dependence of the 3PCF seen 
in simulations which include contributions from realistic aspherical halos with substructure. The left panel shows the result for isosceles 
triangles with r = q = 3' against varying the interior angle tp (see Figure 4 for the triangle geometry). The right panel is for more elongated 
triangles with r = 2q ~ 3' . Here, the scale r ~ 3' is chosen based on the fact that it is in the non-linear regime and unlikely to be affected 
by the resolution of the simulations. All the plots display excellent agreement between the halo model predictions and the simulations. 

4.6 The three-point correlation functions of the shear fields 

We turn to the shear 3PCFs, which are more complex but are important because they are easier to measure from survey data than the 
convergence 3PCF. We will focus on the shear 3PCFs rather than the reduced 3PCF, since there is an ambiguity in defining them for the 

* Note that, even if we employ the halo boundary of rigo as discussed in Figure A4, the bump feature for Qk remains at the same level 
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Figure 8. The reduced 3PCF of the convergence field, Q^, for the ACDM model against triangle configurations parameterized by r, q and iji (see Figure 
4). The left panel shows the result for isosceles triangles with r = g = 3' fixed and varying iji, while the right panel for more elongated triangles with 
r = q/2 = 3'. 
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Figure 9. The shear 3PCFs for equilateral triangles against side length r for the ACDM (upper panel) and SCDM (lower panel) models, as in the bottom 
panel of Figure 6. From symmetry considerations, the four parity-odd functions defined by equation (32) vanish for an i?-mode, and three of the parity even 
functions are equal. Hence, only the results for C+++ and C+ x x shown. Note that the absolute value of x x is plotted, since it becomes negative at 
large scales. The two solid curves are the halo model predictions, using only the 1-halo term. 



shear. The ambiguity arises from the fact that the +/x components are defined with respect to the triangle center (see equation (31) to form 
the 3PCF, hence there is no clear choice for defining the 2PCFs that enter in the denominator of Q. One way to do this is to use or ^.y, x 
defined with respect to the side of the triangle connecting the two vertices, but we will not pursue this. Figure 9 compares the halo model 
predictions with simulation results for the shear 3PCFs for equilateral triangles, as in the bottom panel of Figure 6. For equilateral triangles 
there are only two independent, non-zero 3PCFs: (,+++ and C+x x = Cx + x = Cx x + > while Cx x x = Cx++ = C+x+ = C-t-+x = (see 
§3 and also SL03 and TJ03a). The figure thus shows only the results for C,+++ snd C+x x • The upper and lower panels show the results for 
the ACDM and SCDM models, as indicated. Note that the plots are on a logarithmic scale and the absolute value of ^+ x x is plotted, since 
it becomes negative on large scales. The halo model predictions include only the 1-halo contribution. The (,+++ component carries most of 
the information of the lensing signal for equilateral triangles (ZS03; TJ03a). 

Figure 9 shows that the halo model agrees with the simulation results for angular scales r ^2' and ^ 1' for the ACDM and SCDM 
models, respectively. Whether or not the discrepancy at the smaller scales is genuine is unclear due to the resolution limit of the simulations 
(see §4.2). We found that the shear correlations are more sensitive to the angular resolution of the simulations than the convergence field. 
The agreement extends to scales ^ 10', though the PT contribution to the shear 3PCFs is not included, in contrast to the convergence 3PCF 



© 0000 RAS, MNRAS 000, 1-31 



16 M. Takada & B. Jain 



10" 



8x10"' - 



6x10" 



4x10" 



ojr 2x10" 



-2x10" 



8x10" 



6x10" 



4x10" 



^^2x10" 



-2x10 



parity— even. 



r=q=3 arcmin 3x10 



2x10 




parity— odd 




4x10" 



-2x10" 



4x10 



3x10 



2x10" 



10" 




0.4 0.6 





1 


-10 


-7 1 


2x10 


-7 


3x10 


-7 




0.8 



0.2 



0.4 0.6 



0.8 



Figure 10. The eight shear 3PCFs for the ACDM model against triangle configurations, as in Figure 8. The upper and lower plots show the results for the 
parity-even and -odd functions, respectively. Note that range on the y-axis for the right panel is about two times smaller than in the left panel. The solid curves 
show the halo model predictions for the eight shear 3PCFs, while the symbols are the simulation results as indicated. 



where the PT contribution is dominant at 5'. This agreement is somewhat surprising, but it might be explained as follows. The spin-2 field 
properties of the shear fields cause cancellations between the shear 3PCFs to some extent, which explains why the shear 3PCF amplitude is 
smaller than the convergence 3PCF by an order of magnitude (see Figure 8). The shear 3PCFs appear to be dominated by the coherent shear 
pattern around a halo rather than filamentary structures in the quasi non-linear regime, which are described by PT (at least for triangles that 
are close to equilateral). 

Both the halo model and the simulations show a flattening and decline in at small scales. The halo model predicts the decline at 

slightly smaller scales r ^ 0'.5 compared to the simulations. The origin of this feature may be explained as follows. The shear 3PCFs vanish 
as one goes to zero triangle size (as the three vertices approach the same point), since in this limit it is equivalent to the shear skewness which 
vanishes from statistical symmetry. This limiting behavior is another way to see why the amplitude of the shear 3PCF is smaller than that of 
the convergence. Nevertheless, it is possible that the flattening scale of the shear 3PCF reflects a scale related to the size of typical halos that 
provide the dominant contribution. It will be interesting to study this feature more precisely with higher resolution simulations. 

For C+x X, the agreement between the halo model prediction and the simulation result is not good as for However, the accuracy 

of the simulation results is likely to be worse because of the smaller amplitude of (^+ x x ■ 

The accuracy of the halo model for the shear 3PCFs is tested in greater detail in Figure 10, which compares model predictions with 
simulations for varying triangle configurations for the ACDM model. The upper and lower panels show the parity-even and -odd functions, 
respectively. As can be seen from the lower left panel, both the halo model predictions and the simulation results verify that two parity- 
odd functions vanish: i^xxx ~ Cx++ = for isosceles triangles (see §3). However, the other two parity-odd functions do carry lensing 
information as pointed out by SL03 (also see TJ03a), as the x component is at the vertex bounded by unequal sides. For ?/> = 7r/3 the 
triangle is equilateral and these two functions also vanish. 

The right panel of Figure 10 shows that all the eight functions are non-zero for general triangle configurations (SL03; TJ03a). In contrast 
to the convergence 3PCF shown in Figure 8, the shear 3PCFs display complex configuration dependences and change sign with varying tp. 
These features reflect the detailed structure of the underlying mass distribution as well as properties of the spin-2 field generated by an 
E-mode signal. C+++ peaks around -0 = 7r/3, where the triangle configuration is close to equilateral. For the general triangles shown in 
the right panel, all the eight functions have roughly comparable amplitude. Comparing the left and right panels shows that more elongated 
triangles lead to smaller amplitudes of the shear 3PCFs. These results are to be contrasted with the expectation that the 3D mass 3PCF has 
higher amplitude for elongated triangles on large scales ( ^ lOMpc), reflecting the dominance of anisotropic structures in the perturbative 
regime (Scoccimarro & Frieman 1999; Bemardeau et al. 2002b; also see Figure 5 in TJ03b). The range tp = [0, tt] is shown in our figures; 
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Figure 11. As in the previous figure, for the SCDM model. 



one should keep in mind the relation C(^) = ±C(27r — V') for the lensing i5-mode, where + or — sign is taken for the parity-even or -odd 
functions, respectively (§3; also see Figure 3 in TJ03a ). 

To further check the accuracy of the halo model for different cosmological models. Figures 11 and 12 show the results for the SCDM 
and rCDM models, respectively. Note that the cosmological parameters for the rCDM model are the same for the SCDM model except for 
the shape parameter F = 0.21 (F = 0.5 for the SCDM model). The angular resolution of the two models is the same (Jain et al. 2000). The 
halo model predictions again match the simulation results for all eight functions. A comparison of Figure 1 1 and 12 shows that the amplitude 
of the shear 3PCF depends on the shape of the input linear mass power spectrum. In the context of the halo model picture, this is captured 
by the dependence of the halo mass function on the power spectrum shape on the angular scales we considered. In summary, from Figures 
10-12, one can see that the shear 3PCF amplitudes are sensitive to the cosmological models, but the oscillatory shape of each function is 
quite similar for the three CDM models. 



4.7 Summary: the halo model accuracy for predicting the lensing statistics 

Here we summarize the results shown in Figures 5-12, which investigated the 2PCFs and the 3PCFs of the convergence and shear fields. It 
has been shown that the halo model predictions are in remarkable agreement with the simulation results for all these statistical quantities over 
the angular scales we have considered. In particular, the halo model reproduces the amplitudes and the complex configuration dependences 
for the eight shear 3PCFs. We have used only the 1-halo term and not adjusted any model parameters to get this agreement. This implies 
that the shear 3PCFs result from correlations between the tangential shear pattern around a single NFW profile, as pointed out by TJ03b and 
ZS03. The agreement is striking, since the halo model rests on simplified assumptions of smooth, spherical halos while the halos in CDM 
simulations are aspherical and contain substructure (e.g., Jing Sl Suto 2002). The projection of halos oriented in different directions and the 
nonlocal properties of the shear appear to dilute the effect of some of the detailed structure of halos and make the spherically symmetric 
profile a good approximation in the statistical sense. 

It is not clear whether the agreement remains on smaller angular scales ( ^ 1'), due to the resolution limitation of the simulations we 
have used. It is of great interest to address this issue using higher resolution simulations. On these scales new ingredients may be needed 
for the halo model as well, such as the inclusion of substructure (Sheth & Jain 2002). The agreement we have shown holds for different 
cosmological models (ACDM , SCDM and rCDM models). This implies that the cosmological model dependences can be captured through 
the spherical collapse model and the mass function used in the halo model. These results lead us to conclude that the halo model provides an 
analytical method for predicting higher order lensing statistics with sufficient accuracy for our purposes. 
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Figure 12. As in the previous figure, for tlie rCDM model. The cosmological parameters for the rCDM model are the same as for the SCDM model, except 
for the shape parameter F which is 0.21, while it is 0.5 for SCDM. 



5 MEASUREMENT OF HALO PROFILE PARAMETERS FROM SHEAR CORRELATIONS 

In the following we address the question: how can measurements of the lensing 2PCF and 3PCFs be used to constrain halo model parameters? 
In particular, we focus on the halo profile parameters: the inner slope a for the generalized NFW profile in equation (10) and the halo 
concentration c in equation (12). We will show that forthcoming lensing surveys can put stringent constraints on these parameters. For this 
study we will use the convergence 2PCF and 3PCF for simplicity, although the shear correlations are the direct observable. To use the shear 
3PCFs, it is necessary to combine all the eight 3PCFs. Since the lensing information obtained combining the eight shear 3PCFs are related 
to the convergence 3PCF, the results we obtain should hold as a first approximation. We also note that the statistics of the convergence field 
may be directly measured from future space based surveys such as the one proposed for the SNAP satellite (Jain 2002). 



5.1 Sensitivity of the lensing 2PCF and 3PCF to profile parameters 

In TJ03b, we showed that the 2PCF and 3PCF of the 3D mass distribution at small scales are sensitive to halo profile parameters (see Figures 
12-14 in TJ03b). It was shown that the 3PCF is more sensitive to the halo profile than the 2PCF. This is expected to hold for lensing statistics 
also, since the lensing fields are projections of the 3D mass distribution. 

We derive analytical expressions in Appendix B for the halo convergence for a = 0, 1 and 2 in the generalized NFW profile of equation 
(10). To compute the convergence 2PCF and 3PCF for general a (0 < a < 2) we simply interpolate from the 1-halo term predictions for 
a = 0, 1 and 2, using formulae similar to equation (42) in TJ03b; the interpolation is expected to work to better than 10%. 

Figure 13 shows the sensitivity of the convergence 2PCF (left panel) and 3PCF (right panel) to the inner slope parameter a and to the 
concentration parameter co; in our parameterization c = co(l + z)^^ {M /Mi.)~^ with (3 — 0.13. Increasing a or co steepens the 2PCF and 
3PCF for scales 6 ^ 3', since it increases the density profile in the inner region, r < rvir/c. The curves coincide with each other at large 
scales, because the outer region has the slope for all a. The 3PCF is more sensitive to modifications of the halo profile than the 2PCF; 
this is because of the extra power of the halo profile in the 1-halo term. 

Figure 14 shows the mass range of halos that contribute to the lensing statistics on scales of our interest. The figure shows the depen- 
dences on the maximum mass cutoff in the halo model calculation on the convergence 2PCF and 3PCF. Massive halos with M > IO^Mq 
provide more than 80% of the contribution over the scales we have considered. At smaller angular scales, less massive halos are more 
relevant. One can also see that the 3PCF is more sensitive to massive halos than the 2PCF. 
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Figure 13. The dependences of the 2PCF (left panel) and 3PCF (right panel) of the convergence field on halo profile parameters. The upper and lower dashed 
curves show the halo model predictions with a = 2 and for the inner slope parameter of the generalized NFW profile of equation (10). The upper and lower 
dot-dashed curves are the results for the concentration parameter cq = 15 and 3. The solid curve is the halo model result for our reference model, the NFW 
profile with a = 1, cq = 9. 




Figure 14. The dependences of the 1-halo term contribution to the convergence 2PCF (left panel) and 3PCF (right panel) on the maximum mass cutoff used 
in the calculation. From top to bottom, the five solid curves are the results for a maximum mass of 10^^ , 10^^, 10^'', lO'^^ and 10^^ Mq, as indicated. The 
dashed curve shows the total halo model prediction. Most contributions arise from halos of M > 10^^ Mq on the scales we have considered, and the 3PCF 
is more sensitive to massive halos than the 2PCF. 
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5.2 Covariance of the 2PCF and 3PCF estimators 

To rigorously extract parameter information from cosmic sfiear measurements, it is crucial to consider the covariance between the lensing 
statistics in different bins. This issue for the shear 2PCFs has been investigated by Schneider et al. (2002b). It was shown that there are 
strong correlations between the shear 2PCFs in different bins on small scales. Extending their method, we present analytic expressions for 
the covariance of the convergence 3PCF and 2PCF. Note that the method described below can be extended to compute the covariance for the 
shear 3PCFs, if one accounts for the spin-2 phase factors of the shear fields and the projection operators. 

Following Schneider et al. (2002b), an estimator of the convergence 3PCF from realistic data can be expressed as 

CT(r,q,ii>) = -i;^y^KiK,jKkKjk{r,q,-^), Nuip = y^Aijk{r,q,ip), (38) 

Jvtrip ^ ' ^ ' 

ijk ijk 

where index i denotes source galaxies and A/trip is the number of triplets of galaxies that form a given triangle configuration within the bin 
width. The function Aj^fc (r, g, is the selection function for triplets: it is unity if the three points with indices i, j and k are in the triangle 
configuration and zero otherwise. In the weak lensing limit, the observed convergence field can be expressed as the sum of the lensing 
signal and the noise contamination due to intrinsic elUpticities: k = /t'™'* + Ue. The expectation value of the estimator above is obtained by 
averaging over source ellipticities and performing an ensemble average of the convergence field, denoted by (•••). We will use the notation 

The covariance of the convergence 3PCF is defined as 

cov[a, c;] = ((cr* - a)(cr' - c;)> = (cr*cr*'> - c.C- (m 

where we have simplified the notation so that and Ck denote Cn('", 1, tp) and C«(^', q\ V'')' respectively. Similarly as done in Schneider et 
al. (2002b), the first term on the r.h.s. of the above equation can be rewritten as 

(CrCr> = jv,Hp(r,g,V)^tHp(r',g',^0 E A,,.(r-,g, V')A„.c(r',g',V')(«i«i«/c«a/..Kc). (40) 

ijkabc 

The covariance estimate thus requires an evaluation of the 6-point correlation function of the observed convergence field. If the noise field 
and the convergence field are statistically uncorrelated, the 6-point correlation function can be expressed as 

(KiKjKkKaKbKc) w (Te {SiaSjbSkc + 5 perm.) + [^K{0ij)^i^{0ka)^K.{0bc) + 14 perm.] , (41) 

where ac is the rms of the intrinsic ellipticities " and we have ignored terms of 0{af£^K.) and 0(0-^^^), which are relevant only in the transition 
regime between those in which the first or second term on the r.h.s of equation (41) dominate. In addition, the second term ignores the non- 
Gaussian contribution that is due to the connected parts of the three-, four- and six-point correlation functions, and thus underestimates the 
sample variance. For a more conservative estimate of the non-Gaussian contribution, one might replace the coimected parts of the three-, 
four- and six-point functions with their unconnected parts, providing {nf ■ ■ ■ k^^) = 8 [^Ki6ij)^K{0ka)^Ki0bc) + 14 perms.]. However, 
in this paper we use the above equation for simplicity. This is likely to be a good approximation for two reasons. One, our estimates are 
consistent with the results from simulations which contain the full non-Gaussian contribution, as shown in Figure 15. Moreover the shot 
noise contribution dominates the covariance on sub-arcminute scales, which provide the main constraints on halo profiles. 
Inserting equation (41) into equation (40) and performing an ensemble average over source galaxy positions yields 

Cov[C(r, q, C(r', q' , V-')] = [AC^(r, q, V')]' S{{r, q, ^} - {r' , q' , ^'}) + 7^(r, q, ^, r',q', ^'), (42) 

with 



A<^N{r,q,ip) 



•\/ A'trip 

( ^ 
X — 

Vl' 



lO^degV VlO^ai 

/Alnr\-i/2 / g\-i /AlnqX-i/^ /AtA/ttV'^^ 



V 0.1 J VlV V 0.1 / \ 0.2 



(43) 



(44) 



TZ{r,q,tp,r',q,ip') = J sds j difir J dc^^ [$(r)$(s - r - q)$(qr') + 14 perms.] (s,(/?r,¥'r;r,g,V,r',g',V'). 

where Og and rig are the survey area and the number density of source galaxies, respectively, and Ar, Aq and Atp denote the bin widths 
for the three parameters (r, q, ip) that specify the triangle configuration. The notation used in equation (44) is Smax = -^/Og/Tr, r = 
r(cos <fir, sin (pr), Q = q{— cos(')/' — <Pr), sin(')/) — (pr)) and so on. We have assumed that the survey geometry does not affect covariance 
estimation - a good approximation as long as we consider sufficiently small scales compared to the survey size. The first term on the r.h.s of 
equation (42) denotes the shot noise contribution due to the intrinsic ellipticities, where the function 5{{r, q, 4>} — {r' , q' , V"'}) is defined to 
be unity if r = r',q = q' and ip = ip' within the bin widths, and zero otherwise. The derivation of this term requires an estimate of the triplet 



Van Waerbeke (2000) discussed a more accurate description of the noise field for the convergence, taking into account the smoothing kernel used for the 
reconstruction. 
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Figure 15. Error estimates for measurements of the convergence 2PCF (left panel) and 3PCF (right panel) for the ACDM model. The error estimate includes 
both the shot noise due to intrinsic ellipticities, with rms cr^ = 0.4, and sample variance. Two cases for survey area and number density of source galaxies 
rig are shown. For (Qs.ng) = (1000 deg^, 100 arcmin ), the dotted and dot-dashed curves show the shot noise contribution and the sample variance 
separately: shot-noise dominates on sub-arcminute scales. The solid curve shows the halo model prediction for the 2PCF or the 3PCF. The triangular symbols 
show the sample variance estimated from simulations for f2s = 1000. 



number for a given triangle configuration, which we estimate as A^trip = ngQ^ x ng-Kr{Ar) x ngq{Aq){Atp), where the first, second and 
third factors denote the number of galaxies at the vertices, 1, 2 and 3, respectively (see §4.3). The second term TZ in equation (42) denotes 
the sample variance. Several interesting points made by Schneider et al. (2002b) for the 2PCF hold for the convergence 3PCF as well. (1) All 
the terms are proportional to Q^^' ^^'^ therefore the relative contribution of the terms is independent of survey area, if the area is sufficiently 
large. (2) The sample variance TZ does not depend on the survey particulars such as rig and a^. Further, TZ is independent of the bin widths. 
This implies that combining the 3PCFs in different bins cannot reduce the sample variance. This is indeed verified by ray-tracing simulations. 
(3) The off-diagonal components of the covariance arise only from the sample variance TZ. 

To compute the sample variance TZ for a given cosmological model using equation (42), we first make a table of the model predictions 
of the convergence 2PCF as a function of the separation angle. Then, we can use the same table to compute the sample variance between the 
two convergence 3PCFs of any triangle configurations. Therefore, this method is much more tractable than a direct implementation including 
contributions from the three-, four- and six-point functions, where we have to account for their configuration dependences in the integration 
of equation (44). An alternative way to estimate the covariance is to use ray- tracing simulations. However, to do this rigorously requires an 
adequate number of independent realizations, since sample variance is large for the higher-order moments. The PTHalos method recently 
proposed by Scoccimarro & Sheth (2002) can be a powerful tool for such an approach. 

Likewise, we can derive the covariance of the convergence 2PCF as 

Cov[C(r),C(/)] = [AU(r)]^5(r-r') +7^2pt(r,/) (45) 
with 

7^2pt(r,r') = sdsj d^rj^ d^:.[as)a\s + r' -r\)+a\s + r'ms-r\)], (47) 

where the number of pairs are estimated as A'^pair = UgQs x ngTTr{Ar). 

Figure 15 plots the square root of the diagonal component of the covariance matrix for the convergence 2PCF (left panel) and 3PCF 
(right panel) for the ACDM model. This is an estimate of the error on the 2PCF and 3PCF measurements from a lensing survey. We consider 
two cases, specified by the survey area fls and the number density of source galaxies Ug. The dashed and broken curves show the results 
for {Qs, Ug) = (1000, 100) and (200, 40) in units of degree^ and arcmin"'^, respectively. The former is expected from future imaging 
survey, while the latter applies for a survey like the CFHT Legacy Survey, which has just begun. For (ils, Ug) = (1000, 100), the dotted and 
dot-dashed curves show the shot noise contamination and the sample variance separately, implying that the shot noise provides the dominant 
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contribution at small scales 1'. The triangle symbols denotes the simulation results for the sample variance for Os = 1000 degree^. Note 
that the simulation result is computed from 36 realizations of a 11.7 degree^ simulated map, and is then scaled as oc ^ts ^^'^ . Out analytic 
estimates are consistent with the simulation results, within the rather large error bars of the latter (not plotted to preserve clarity). 

The solid curve in each panel of Figure 15 shows the halo model prediction for the 2PCF or the 3PCF. Comparing the solid curves with 
the error estimate gives the signal-to-noise (S/N) ratio for measuring the 2PCF and 3PCF. Clearly these survey parameters will allow for 
measurements of shear correlations at high significance, even on sub-arcminute scales. We note that the S/N estimate for the 3PCF is only 
for one triangle configuration - we can combine the 3PCFs from different configurations to improve the S/N at these scales. 



5.3 Constraints on a and co 

Next we apply the covariances computed above to demonstrate how combined measurements of the 2PCF and 3PCF can be used to constrain 
parameters of the halo profile. We use the standard statistic, expressed for our case as: 

= X2pt + X3pt' (48) 
where 

X2pt = ^(|.-^i)[Cov(0]^-'fe-6), 

i<j 

X3pt = ^{Ci-Ci)[GoY{0]7/{Q-Q), (49) 

where ^ and C. denote the 2PCF and 3PCF for the fiducial model and [Cov]^-^ denotes the inverse of the covariance matrix. The index i runs 
among different bins for ^ and ^. In equation (48), we have ignored the covariance between the 2PCF and 3PCF for simplicity. We make this 
approximation because only the sample variance contributes to this covariance (the shot noise contribution vanishes as it involves the fifth 
power of the intrinsic ellipticities), which is small on sub-arcminute scales as discussed above. 

So far we have used the parameters r, q and i/> to describe triangle configurations. To perform the fitting we employ an alternative 
set of parameters used in the literature (e.g., Peebles 1980): 

r = 0i2, ■"=— , v= — , (50) 

fl2 fl2 

with the condition 9i2 < 623 < 631, which imposes the constraints w > 1 and < v < 1. Different sets of r, u and v correspond to different 
triangles, so we do not have to worry about double-counting. 

We treat the inner slope parameter a and the halo concentration normalization co as free parameters. We fix /3 = 0.13 (the mass 
dependence of the concentration) and we set the other model parameters to be those for the ACDM model. We discuss below why this choice 
is a good first attempt at the use of correlation statistics to measure halo profiles. As shown in Figure 13, these profile parameters are sensitive 
to the lensing statistics at small scales ^ 2'. We therefore restrict ourselves to these scales. We consider 10 logarithmic bins in r = [0'.12, 2'] 
with the bin width Ar/r — 0.1. For the 3PCF, we consider 5 bins each for u and v. u = 1, 2, 3, 4, 5 and w = 0.1, 0.3, 0.5, 0.7, 0.9. Thus 
we use 10 bins for the 2PCFs and 250 for the 3PCF. How the binning affects the fitting of model parameters must be carefully examined, 
to avoid over- or imder-estimating the constraints (e.g, Scoccimarro & Frieman 1999). Since we correctly take into account the covariance 
for the 2PCF and 3PCF, including the off-diagonal components, we avoid over-constraining the parameters. The triplet number used for the 
shot noise evaluation is estimated for the binning we consider as A^iip = "gfi„ x ivrAr x qrAuAv with Au = 1 and A;; = 0.2. To save 
computational expense for the sample variance of the 3PCF, we ignore the dependence on the v parameter - we compute the sample variance 
for different bins of r and u, but with v = 0.5 fixed, resulting in 50 x 50 computations for the sample variance. This is adequate, since the 
configuration dependence is weak as shown in Figure 8. 

Figure 16 shows contour plots of the constraints on the inner slope parameter a and the halo concentration cq. We consider S7s — 200 
degree^ and ng — 40 arcmin^^ for the survey area and the number density of source galaxies, respectively. The left panel shows the 
constraints if we use the measurement of the 2PCF only, while the right panel shows the constraints from combined measurements of the 2PCF 
and the 3PCF. The cross symbol denotes our fiducial model of (a, co) = (1, 9), which is the NFW profile with the concentration given by 
Bullock et al. (2001). For a Gaussian probability distribution function for the 2PCF and 3PCF, the three shaded regions correspond to Ax^ ~ 
2.30, 6.17 and 11.8 corresponding to 68.3%, 95.4% and 99.73% confidence levels (C.L.), respectively. One can see that the constraints 
on a and co are degenerate: the effect of increasing (decreasing) a on the 2PCF and 3PCF is compensated by decreasing (increasing) 
Co. Nevertheless, the right panel shows that the 3PCF measurement can tighten the constraints, implying that the 3PCF provides additional 
information on these parameters, although the parameter degeneracy is not broken. As a result, this type of survey can be used to put stringent 
constraints on a and co within 5% level (95% C.L.), if one of the parameters is fixed, and systematic uncertainties in the measurements and 
the theoretical model do not dominate. 

Figure 17 shows the result expected from a different type of survey from the one in the previous figure: a deep, small area survey, with 
(fis, ng) = (20, 100). For this case, the 3PCF can substantially tighten the constraints provided from the 2PCF. However, the degeneracy 
between a and co remains. Recently, Takada & Hamana (2003) proposed that a joint measurement of the magnification statistics and the 
cosmic shear could be used to break the degeneracy or at least put upper bounds on a and co, since the amplitude of the magnification 
correlation is more sensitive to an increase of a and co than the cosmic shear correlation. This arises from the non-linear relation between 
the magnification and shear, given by// = \{1 — k)"^ — ■y^\~^. 
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Figure 16. Contour plots of the constraints on the inner slope parameter a and the halo concentration cq expected from a survey with area Hs = 200 
degree^ and Ug = 40 arcmin"^. The left panel shows the result if we use only the 2PCF measurement, while the right panel shows the result from 

combined measurements of the 2PCF and 3PCF. The cross symbol denotes input fiducial model of (a, co) = (1.0, 9.0). To perform the fitting, we have 
used measurements on scales [0'.2, 2'] and combined 250 triangle configurations for the 3PCF (see text for details). The three shaded regions represent 
Ax^ = 2.30(68.3%), 6.17(95.4%) and 11.8(99.7%), respectively. The 3PCF measurement can tighten the constraints, reflecting the fact that it carries 
lensing information complementary to the 2PCF. 



6 DISCUSSION 

In this paper, we have developed the halo model for computing the higher order correlations of the cosmic shear field, extending the real-space 
dark matter halo approach developed in TJ03b. A detailed investigation of the three-point correlation functions (3PCF) of the convergence and 
shear fields with respect to the size and configuration dependence of triangles has been presented. Our method provides an accurate, anal3ftical 
way of computing the 3PCFs with little computational expense. We have focused on the eight shear 3PCFs defined from combination of the 
+/ X projections of the shear fields at each vertex of a given triangle (SL03, ZS03 and TJ03a). The shear 3PCFs defined in this way have 
characteristic properties that help disentangle the lensing i?-mode from the the B-mode due to possible systematic errors and other non-linear 
effects (see §3). They are a direct probe of the gravitational clustering of the mass distribution and can provide an independent test of the 
CDM paradigm of structure formation. 

We have carefully checked the accuracy of our model by comparing the predictions with ray-tracing simulation results. We paid par- 
ticular attention to the triangle configuration dependences of the 3PCFs, since our halo model uses simple spherically symmetric profiles. 
We find excellent agreement over the angular scales and models we have considered, as shown in Figures 5-12. The halo model reproduces 
the complex configuration dependences for the eight shear 3PCFs, as well as their amplitudes. The agreement is found for plausible model 
ingredients: mass function (Sheth & Tormen 1999) and NFW halo profiles with recent prescriptions for the halo concentration (Bullock et al. 
2001). We chose the best parameter values identified in the literature and did not adjust any parameters. 

On scales ^ 1' the 3PCF can be used to break degeneracies in cosmological parameters, in particular in the CTs-fimo determinations so 
far made from the 2PCF measurement. Figures 6-12 show a clear dependence of the lensing 3PCFs on the cosmological models. In addition, 
on large scales the 3PCF can constrain primordial non-Gaussianity, which can be separated from the gravitationally induced signal using its 
dependence on scale (R. Scoccimarro, private communication). In practice, measuring the shear 3PCFs is more feasible than the convergence 
3PCF, since obtaining the convergence requires a non-local reconstruction from the observed shear field. Detections of shear three-point 
moments have recently been reported (Bernardeau et al. 2002a; Pen et al. 2003). Thus, current survey data are likely to already allow for 
shear 3PCF measurements (see Figure 15 and Figure 4 in TJ03a for theoretical justification). Our method provides the only well-tested 
analytical approach to interpret the measured signals in terms of cosmological parameters. 

A second application we have proposed is to use shear statistics on sub-arcminute scales to constrain halo profile properties. Forthcoming 
lensing surveys promise to measure the sub-arcminute signals with high significance (see Figure 15; also see Figures 17 and 18 in Van 
Waerbeke et al. 2001b for a measurement of the shear 2PCFs). This will open a new window in the use of shear correlation functions, beyond 
the determination of cosmological parameters. The n-point correlation functions on sub-arcminute scales arise mainly from correlations 
between the lensing fields around a single halo of M > lO^^M© (see Figure 14). The halo model allows us to interpret the measured 
signals in terms of the halo profile properties. The inner slope of the generalized NFW profile and the halo concentration are sensitive to the 
amplitudes of the lensing 2PCF and 3PCF on these scales (see Figure 13). 

We have demonstrated how combined measurements of the 2PCF and 3PCF can put stringent constraints on halo profile properties 
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Figure 17. As in the previous figure, but for a survey with Clg = 20 degree^ and ng = 100 arcmin ^ . The results show the value of increased number density 
of galaxies that can be achieved with a deeper survey. For such a survey the 3PCF provides a much larger improvement in the constraints. 



from forthcoming lensing data. This was done by taking into account the covariance for the 2PCF and 3PCF measurements, which contains 
contributions from the shot noise due to the intrinsic ellipticities and sample variance. For example. Figure 16 shows that a survey with 
parameters similar to those of the CFHT Legacy survey can constrain the iimer slope parameter with 5% accuracy (95% C.L.) if the halo 
concentration is fixed. Figure 17 shows the dramatic improvements possible with a deeper survey; on the basis of these figures, it follows 
that a deep survey with area ^ 200 square degrees and source number density ~ 100 per square arcminutes can achieve an accuracy of 1% 
in density profile parameters. The use of the 3PCF is critical in being able to achieve this accuracy, and the use of the four-point function 
would also be valuable if it could be measured. In this paper we have ignored the effect of uncertainties in the mass function and its possible 
degeneracy with the inner slope and concentration. White (2002) suggested adjustments to the parameters in the Sheth-Tormen mass function 
from fitting to simulations. We confirmed that this modification alters the halo model predictions for lensing correlations over angular scales 
( ^ 5') where the 1-halo term is relevant. Therefore, measuring the lensing 2PCF and 3PCF can be similarly used to constrain the shape of 
halo mass function over a mass range of 10^^ — IO^^Mq. This would be complementary to other methods such as cluster counts. Breaking 
the degeneracy in the halo model parameters will require either using the different 3PCFs of the shear, which we have not done here, or input 
from other methods such as the convergence reconstruction approach discussed below, or a statistical study of strong lensing. X-ray and SZ 
effect on cluster/group scales. These are important issues to be addressed. 

A fundamental result from CDM simulations is that the density profiles of halos are universal across a wide range of mass scales (e.g., 
NFW). Therefore, applying our method to different length scales would sample halo profiles on different mass scales and offer a powerful test 
of the CDM paradigm. The halo model formalism can also be extended to include the effects of substructure, currently a subject of study due 
to a possible conflict between CDM theory and observations (e.g. Sheth & Jain 2002; Dalai & Kochanek 2002). Substructure and triaxiality 
of halos, discussed further below, would increase the amplitude of correlation functions on the smallest scales. Since the two-, three-, and 
four-point functions scale differently with these effects, a careful study is merited to develop the correlation function approach as a probe 
of different small scale effects. Using lensing correlations on scales of 0.1 — 2 arcminutes, halos with masses of 10^^ — M^^Mq will be 
probed. 

The method of constraining halo profile properties from shear correlations is complementary to the approaches of Reblinsky & Bartel- 
mann (1999), Dahle, Hannested & Sommer-Larsen (2003), White et al. (2002), Miyazaki et al. (2002) and Padmanabhan et al. (2003). They 
rely on the reconstruction of the convergence field from observed ellipticities (see Dahle et al. 2003 and Miyazaki et al. 2002 for the imple- 
mentation to actual data). Halo profiles are constrained by looking at convergence profiles around individual peaks in the convergence map. 
However, there are some limitations to this method. The first is lensing projection effects - we cannot accurately measure properties of the 
primary halo due to superposition with a void region or less massive halos along the same line of sight, even if the redshift of the lens halo is 
available (White et al. 2002; Padmanabhan et al. 2003). The other limitation is the angular resolution of the convergence reconstruction. In 
practice, reconstructing the convergence field on a given patch of the sky requires an averaging of the observed ellipticities over an adequate 
nimiber of source galaxies to reduce the noise contamination as well as enhance the contrast of the signals arising from halos. For plausible 
survey parameters the reconstruction resolution is larger than an arcminute. The reconstructed convergence around a halo is thus smoothed, 
which makes it difficult to see the inner region in the halo profile. Padmanabhan et al. (2003) concluded that the inner slope of NFW profiles 
caimot be constrained using convergence maps. 

The resolution limitation of convergence maps can be offset by stacking clusters and by follow-up deeper weak lensing observations 
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and other methods such as the SZ effect, X-ray observations and strong lensing. It is not an easy task however, and is subject to biases 
associated with identifying centers and mass scales with peaks. The strength of the correlation function method we have discussed is that it 
allows for constraints on the inner halo profile statistically, even though individual halos are not resolved on these scales. Our approach treats 
the data objectively, while taking a parameterized approach to the modeling. This appears a sensible approach at a time when cosmological 
parameters are well constrained and we have a broad understanding of halo properties. 

The real-space halo model formulation, developed in TJ03b and this paper, does not rely on a model of the 3D non-linear power 
spectrum. This fact leads to interesting applications for lensing statistics. We can directly compute any n-point correlation functions of 
lensing observables, such as the reduced shear field 5 = 7/(1 — /t) and the lensing magnification n = 1/[(1 — k)^ -I- 7^] (Takada & Hamana 
2003). This can be done merely by replacing 7m in equation (20) with qm or iim for a given halo of mass M, respectively. The resulting 
prediction is exact in the sense that it fully accounts for the non-linear contribution of lensing. So far, the cosmological interpretation of 
lensing two-point statistics has been made using theoretical predictions computed from a model of the 3D power spectrum, which requires 
a perturbative approach, such as setting g Ki -f or n Ki 1 + 2/t (e.g.. Van Waerbeke et al. 2001b; McKay et al. 2001; Guzik & Seljak 2002; 
Benitez & Martinez-Gonzalez 1997; Gaztaiiaga 2003). 

There are some uncertainties we have ignored in the halo model formulation. We have employed a spherically symmetric halo profile, 
while CDM simulations show that halos are triaxial (Jing & Suto 2002). In addition, high-resolution simulations have also shown that about 
10% of mass distribution in a halo is in small sub-clumps (e.g, Ghigna et al. 2000). We have shown that the halo model computed from the 
NFW profile matches the simulation results for all the lensing statistics we have considered. However, it is unclear whether this agreement 
holds on sub-arcminute scales. The 3PCF can be the lowest order statistical quantity to probe the detailed mass distribution of halos via 
its dependence on triangle configurations. Substructure and triaxiality are expected to have different scale and configuration dependences 
than changes to the halo parameters that we have considered. It is an interesting problem for future work to work out in detail the different 
effects that emerge on scales of order ~ 0.1 arcminute. For the applications described above, it is important to test the analytical predictions 
with ray-tracing simulation with higher resolution. This is a pressing need and a challenge for future niunerical work. The PTHalos method 
developed by Scoccimarro & Sheth (2002) could be a powerful alternate tool for such a study. 

We would like to thank T. Hamana for kindly providing his ray-tracing simulation data and for discussions. We also thank R. Scocci- 
marro, R. Sheth, M. Jarvis and G. Bernstein for valuable comments and discussions and thank R Schneider for valuable discussions on the 
parity transformation properties of the shear three-point functions. This work is supported by NASA grants NAG5-10923, NAG5-10924 and 
a Keck foundation grant. 



APPENDIX A: VALIDITY OF THE REAL-SPACE HALO APPROACH FOR SHEAR STATISTICS 

In this appendix, we address the validity of the real-space halo model developed in §2.4. 

Figure Al plots the shear 2PCF, £,-,{9) = (7 ■ 7*), against separation angle 6 for the ACDM model and source redshift Zs = 1. 
The thick dashed and dotted curves are the real-space and Fourier-space halo model predictions for the 1-halo term, which are computed 
using equations (20) and (24), respectively. Note that we have used the NFW profile. There is excellent agreement over the angular scales 
we have considered, verifying that the real-space model formulated in §2 is equivalent to the Fourier-space model. This agreement is quite 
encouraging, since the real-space halo model for the shear field contains an infinite integration range to account for the non-local property of 
the shear fields. To more explicitly demonstrate the importance of the integration range, the thin dashed curve is the real-space halo model 
prediction when the integration range /d^s is confined to the virial region, as is done for the convergence field (see equations (19)). It 
significantly underestimates the 1-halo term. Figure 3 in TJ03b also demonstrates that the real-space halo model for the convergence field is 
equivalent to the Fourier-space halo model. These results lead us to conclude that the real-space halo model formulation for the shear and 
convergence fields has been made in a self-consistent way, in agreement with the Fourier-space model well studied in the literature (e.g., 
Cooray & Sheth 2002). 

The thin dotted curve denotes the 2-halo term, which provides the dominant contribution to the 2PCF for 9 ^ A' . For comparison, the 
thin solid curve is the result for computed from the linear mass power spectrum. It very slightly overestimates the 2-halo term which 
includes the biasing of halos. The 2-halo term thus captures the clustering properties in the linear regime. 

The accuracy of the halo model for predicting the lensing 2PCFs can be seen by comparing the dot-dashed and thick solid curves (the 
detailed comparison with ray-tracing simulation results will be presented below). The dot-dashed curve shows the total halo model prediction 
(1- plus 2-halo terms), while the thick solid curve is the result computed from the fitting formula proposed in Smith et al. (2002) (hereafter 
Smith02). The Smith02 formula was calibrated to reproduce the non-linear mass power spectra from high resolution iV-body simulations 
for various cosmological models. As can be seen, the halo model prediction matches the Smith02 result within 10% at O'.l < < 30'. The 
reliable range of the Smith02 formula, is fc ^ 10 /iMpc^^ in the mass power spectrum. Since the lensing field is the projected field of the 3D 
mass distribution, the valid range roughly corresponds to the angular scale 6 ^ 27r/[A:maxrfA(z = 0.4)] ~ 0!6, where 2 — 0.4 is close to the 
peak redshift of the lensing projection function for the source redshift Zs = 1 (e.g., see Figure 4 in TJ02). Therefore, it is unclear whether or 
not the discrepancy between the halo model and the Smith02 result at ^ ^ O'.l is genuine. 

It is worth mentioning why we use the fiducial model (co, /3) — (9, 0.13) for the halo concentration in this paper, since we used a 
steeper mass slope l3 = 0.2 in TJ02 and TJ03b. The Smith02 formula predicts a steeper k slope of the non-linear power spectrum P{k) than 
predicted from the fitting formula of Peacock & Dodds (1996) (hereafter PD), which has been widely used in the literature. For the halo 
model, the fe-slope is determined by the halo profile parameters, the slope of the mass function, and the input linear power spectrum (Seljak 

1" The PD formula is based in part on the physically motivated stable clustering hypothesis, which allows one to analytically predict the behavior of strongly 
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Figure Al. The two-point correlation function of the shear field, S;'y{0) = {f ■ "7*), vs. separation angle 9. The thick dashed curve shows the real-space halo 
model prediction for the 1-halo term, computed using equation (20), while the thick dotted curve shows the Fourier-space halo model prediction. The sum with 
the 2-halo term, denoted by the thin dotted curve, leads to the halo model prediction for the total power (lh-H2h terms) shown by the dot-dashed curve. For 
compaiison, the thick and thin solid curves show the predictions computed from the Smith02 fitting formula for the non-linear mass power spectrum and the 
linear power spectrum, respectively. The thin dashed curve is the real-space halo model prediction when the integration range is confined to the virial radius. 



2000; Ma & Fiy 2000a,b,c; TJ03b). The halo model with P — 0.13 gives a slightly better fit to the fc-slope of the Smith02 power spectrum 
than with (3 — 0.2. It is reassuring that we find agreement with N-body simulations with plausible ingredients for the halo model, each of 
which is supported by independent studies: the halo profile (NFW), the halo concentration (Bullock et al. 2001) and the halo mass function 
(Sheth & Tormen 1999). These studies explored halo properties by resimulating regions containing halos in a larger scale simulation with a 
higher resolution simulation. Hence, one advantage of the halo model is that it can be easily refined by incorporating results from different 
A'^-body simulations with different sizes. 

Figure A2 compares the theoretical predictions for ^.^,+ (6) (left panel) and (9) (right panel), as done in the previous figure. The 
real-space halo model result for the 1-halo term (thick dashed curve) again agrees with the Fourier-space result (thick dotted curve) for both 
5.y,+ and ^^.x . The agreement is as a result of the integration of the phase factors e^e^ in equation (21) for the real-space case. Comparison 
of the dot-dashed curve with the solid curve in each panel shows that the halo model prediction for the total 2PCF matches the Smith02 
result. Note that, although the 2-halo term (or the linear theory prediction) for ^^..x overestimates the amplitude at large scales ^ 5', the 
agreement of the halo model with the Smith02 result at these scales results from sum of the negative 1-halo term plus the 2-halo term. 

Figure A3 compares theoretical predictions for from various models of non-linear gravitational clustering. The thick solid, dashed 
and dot-dashed curves are the results from the Smith02 formula, the PD formula and the halo model, respectively. The lower panel shows 
the deviation relative to the Smith02 result (thick solid curve) and shows agreement among these models within 10% over the scales we have 
considered. The Smith02 prediction is higher than PD by 0— 10% at 1 < S < 10'. This discrepancy also exists for the aperture mass variance 
(Schneider et al. 1998), which has been used to disentangle the E/ B-modes from the actual measurement (e.g, Van Waerbeke et al. 2001b; 
Jarvis et al. 2003). Most of the cosmic shear analysis to constrain cosmological models has been made by comparing the measured signals 
with the theoretical predictions computed from the PD formula. Hence, the discrepancy shown implies that as might be overestimated by up 
to ~ 5%, if the constraint is obtained from 6 < 10'. On the other hand, the halo model slightly overestimates the Smith02 and PD results on 
6 > 10', where the 2-halo term yields the dominant contribution. This is because the standard implementation of the halo model imposes the 
condition that the 2-haIo term reproduce the shear 2PCF computed from the linear power spectrum on large scales - this cannot reproduce 
the suppression seen in the realistic non-linear power spectrum over the transition scales between the non-linear and linear regimes. Finally, 
the sensitivity of the Smith02 result to the input linear power spectrum is shown by the thin solid curve, which uses the transfer function 
proposed in Eisenstein & Hu (1999), which is more accurate than the BBKS plus Sugiyama model. The comparison shows a difference less 



non-linear clustering. The Smith02 results display a weak violation of the stable clustering hypothesis. This violation can occur for the halo model as well (Ma 
& Fry 2000b, c; also see Scoccimarro et al. 2001; TJ03b). 



© 0000 RAS, MNRAS 000, 1-31 



Three-Point Correlations in Lensing 27 





6 (arcmin) 

Figure A3. Comparison of model predictions for g-y. Tlie solid, dashed and dot-dashed curves are the results from the Smith02 formula, the PD formula and 
the halo model, respectively. The lower panel explicitly shows the deviation of each model prediction from the Smith02 model. Although all the results agree 
with each other within 10%, the Smith02 result for (^-^ is higher than the PD result by - 10% at 0(1 J: 6 10' (see text for details). Note that the Smith02 
result is computed using the BBKS transfer function. To demonstrate the effect of the input linear power spectrum, the thin solid curve is the Smith02 result 
when we use the fitting formula of Eisenstein & Hu (1999), which is more accurate. 
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Figure A4. The figure shows how the halo boundary condition employed in the halo model calculation affects the prediction for . The thick dot-dashed and 
dashed curves are the results when we employ the NFW profile truncated at the virial radius and at the radius rigo, respectively. The radius rigo is defined 
so that mean density enclosed by a sphere with radius rigo is 180 times the background density. The upper dotted curve shows the halo model prediction if 
one uses the expression in Bartelmann (1996) for the shear profile to compute the 1-halo term. For reference, the thin solid curve shows the Smith02 result in 
Figure Al. 



than 3%. This sensitivity also holds for the halo model prediction. Hence, we can safely use the BBKS plus Sugiyama model to compute the 
halo model predictions. 

Figure A4 explores how a modification of the halo boundary condition used in the halo model calculation affects the model prediction. 
As stated below equation (6), the Sheth-Tormen mass function (6) tends to better fit the mass function measured from A'^-body simulations, 
if one employs the halo mass estimator, Migo, enclosed within a region of the overdensity Aigo {— 180), than the virial mass estimator (e.g., 
White 2002). Therefore, one possible modification of the halo model is to employ the halo profile truncated at the radius rigo If we assume 
that the mass distribution in a halo follows an NFW profile up to rigo, we can obtain the relation between M and Migo, which allows us to 
re-express all the relevant quantities in terms of the new mass Migo, as demonstrated in Hu & Kravtsov (2003). In addition, for consistency 
with the simulation results in White (2002), we employ the parameters of a = 0.67 and p — 0.3 for the Sheth-Tormen mass function (6) in 
the halo model calculation (see Table 2 in White 2002). The thick dashed curve plots the halo model prediction for ^.^ for the halo boundary 
ngo, while the thick dot-dashed curve is the result for the virial boundary. As can be seen, the two results are close, although the halo model 
of rigo better matches the Smith02 result denoted by the thin solid curve over a range of the scale 0'.5 ^9^5'. This is mainly due to the 
enhancement of the 1-halo term in the halo model with rigo, since rigo > '"vir (At, = 334 > 180) for the ACDM model, hence the 1-halo 
term covers a larger range than the virial region. The dotted curve shows the halo model result when one employs the expression for the shear 
profile 7a/ in Bartelmann (1996) to compute the 1-halo term. The expression is derived by the line-of-sight projection of the NFW profile, 
allowing it to extend to infinite radius. It substantially overestimates the amplitude for over the scales we have considered. Hence, if one 
intends to account for the mass contribution outside the virial region, it is necessary to first modify the mass defined within the new halo 
boundary and then to modify the halo model ingredients in terms of the new halo mass. This could improve the halo model accuracy over the 
transition scales between the non-linear and linear regimes, since the mass distribution outside the virial region is relevant for the quasi-linear 
regime (5^1; also see the discussions around Figure 8 in TJ03b). We have confirmed that the boundary condition rigo slightly improves 
the agreement between the halo model prediction for the lensing 3PCFs and the simulation results. However, in this paper, we implement the 
virial boundary for simplicity. 

To summarize the results shown in Figures A1-A2, we have shown that our real-space halo model formulation for cosmic shear statistics 
is self-consistent with the Fourier-space halo model. A great advantage of the real-space halo model is that it enables us to analytically 
compute any n-point correlation functions for both the convergence and shear fields on small angular scales, without additional computational 
effort compared to the two-point function. 



Note that, as long as the halo mass function is then considered a function of Migo, mass conservation and the normalization of the mass function are not 
violated 
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APPENDIX B: CONVERGENCE FIELD FOR A GENERALIZED NEW PROFILE 



In this appendix, we present the convergence field around a generalized NFW profile with a 7^ 1 given in equation (10). 

For general a the convergence field cannot be analytically computed. However, if a = and 1, we obtain analytical expressions to 
evaluate km in equation (13) in terms of Em, which is given by 

Y,M{e) = ^^^F{cd/9.ir), (Bl) 



with 
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(x<l) 
(x = l) 

(1< X < c), 



(a;<l) 
(cc = l) 
(1 < a; < c), 



(B2) 



(B3) 



(B4) 



for a — 2, respectively. The factors /o and /2 in the above equations are f^^ = — c(2 + 3c) / (2(1 + cf) + ln(l + c) and /j"^ = ln(l + c). 
We again note that the convergence fields above are defined from the generalized NFW profile truncated at the virial radius, leading to 
Sm(^) = for 6 > 6'vir. These expressions are used to address the sensitivity of the lensing statistics to the inner slope parameter a. 
Similarly, we can derive analytical expressions for the shear profiles for a = and 2, as in equation (16). 
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